vtkNetCDFCAMReader.cxx 27 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkNetCDFCAMReader.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkNetCDFCAMReader.h"

#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkFieldData.h"
#include "vtkFloatArray.h"
21
#include "vtkIncrementalOctreePointLocator.h"
22 23
#include "vtkInformation.h"
#include "vtkInformationVector.h"
24
#include "vtkMath.h"
25
#include "vtkNew.h"
26 27
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
28
#include "vtkPolygon.h"
29 30 31 32
#include "vtkSmartPointer.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnstructuredGrid.h"

33
#include <vector>
34 35
#include <vtk_netcdfcpp.h>

36 37 38 39 40 41
namespace
{
// determine if this is a cell that wraps from 360 to 0 (i.e. if it's
// a cell that wraps from the right side of the domain to the left side)
  bool IsCellInverted(double points[4][3])
  {
42 43 44 45 46
    // We test the normal 3 points at a time. Not all grids are well-behaved
    // i.e. consistenly use 0 or 360. We've had grid where 3 points on the left
    // side, and just 1 on the right. Just checking the first 3 points (which is
    // what ComputeNormal() does, we may (and do) miss a few cells.
    // See BUG #0014897.
47
    double normal[3];
48 49 50 51 52 53
    vtkPolygon::ComputeNormal(3, points[0], normal);
    if(normal[2] > 0)
      {
      return true;
      }
    vtkPolygon::ComputeNormal(3, points[1], normal);
54 55 56 57 58 59
    if(normal[2] > 0)
      {
      return true;
      }
    return false;
  }
60 61 62 63 64 65

  template <class T>
  inline bool IsZero(const T& val)
    {
    return std::abs(val) < std::numeric_limits<T>::epsilon();
    }
66 67
}

68 69 70 71 72 73
vtkStandardNewMacro(vtkNetCDFCAMReader);

//----------------------------------------------------------------------------
vtkNetCDFCAMReader::vtkNetCDFCAMReader()
{
  this->FileName = NULL;
74
  this->CurrentFileName = NULL;
75
  this->ConnectivityFileName = NULL;
76
  this->CurrentConnectivityFileName = NULL;
77 78
  this->PointsFile = NULL;
  this->ConnectivityFile = NULL;
79
  this->SingleLevel = 0;
80 81
  this->TimeSteps = NULL;
  this->NumberOfTimeSteps = 0;
82 83 84 85 86 87 88 89
  this->SetNumberOfInputPorts(0);
  this->SetNumberOfOutputPorts(1);
}

//----------------------------------------------------------------------------
vtkNetCDFCAMReader::~vtkNetCDFCAMReader()
{
  this->SetFileName(NULL);
90
  this->SetCurrentFileName(NULL);
91
  this->SetConnectivityFileName(NULL);
92
  this->SetCurrentConnectivityFileName(NULL);
93 94 95 96 97 98
  delete this->PointsFile;
  this->PointsFile = NULL;
  delete this->ConnectivityFile;
  this->ConnectivityFile = NULL;
  delete []this->TimeSteps;
  this->TimeSteps = NULL;
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
}

//----------------------------------------------------------------------------
int vtkNetCDFCAMReader::CanReadFile(const char* fileName)
{
  NcFile file(fileName, NcFile::ReadOnly);
  return file.is_valid();
}

//----------------------------------------------------------------------------
void vtkNetCDFCAMReader::SetFileName(const char* fileName)
{
  vtkDebugMacro(<<" setting FileName to " << (fileName?fileName:"(null)") );
  if ( this->FileName == NULL && fileName == NULL)
    {
    return;
    }
  if ( this->FileName && fileName && (!strcmp(this->FileName,fileName)))
    {
    return;
    }
120 121 122 123
  delete this->PointsFile;
  this->PointsFile = NULL;
  delete [] this->FileName;
  this->FileName = NULL;
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
  if (fileName)
    {
    size_t n = strlen(fileName) + 1;
    char *cp1 =  new char[n];
    const char *cp2 = (fileName);
    this->FileName = cp1;
    do
      {
      *cp1++ = *cp2++;
      }
    while ( --n );
    }
   else
    {
    this->FileName = NULL;
    }
  this->Modified();
}

//----------------------------------------------------------------------------
void vtkNetCDFCAMReader::SetConnectivityFileName(const char* fileName)
{
  vtkDebugMacro(<<" setting ConnectivityFileName to "
                << (fileName?fileName:"(null)") );
  if ( this->ConnectivityFileName == NULL && fileName == NULL)
    {
    return;
    }
  if ( this->ConnectivityFileName && fileName &&
       (!strcmp(this->ConnectivityFileName,fileName)))
    {
    return;
    }
157 158
  delete this->ConnectivityFile;
  this->ConnectivityFile = NULL;
159
  delete [] this->ConnectivityFileName;
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
  if (fileName)
    {
    size_t n = strlen(fileName) + 1;
    char *cp1 =  new char[n];
    const char *cp2 = (fileName);
    this->ConnectivityFileName = cp1;
    do
      {
      *cp1++ = *cp2++;
      }
    while ( --n );
    }
   else
    {
    this->ConnectivityFileName = NULL;
    }
  this->Modified();
}

179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
//----------------------------------------------------------------------------
int vtkNetCDFCAMReader::RequestInformation(
  vtkInformation* vtkNotUsed(reqInfo),
  vtkInformationVector** vtkNotUsed(inputVector),
  vtkInformationVector* outputVector)
{
  if(this->FileName == NULL)
    {
    vtkWarningMacro("Missing a file name.");
    return 0;
    }

  if(this->CurrentFileName != NULL &&
     strcmp(this->CurrentFileName, this->FileName) != 0)
    {
    delete this->PointsFile;
    this->PointsFile = NULL;
    this->SetCurrentFileName(NULL);
    }
  if(this->PointsFile == NULL)
    {
    this->PointsFile = new NcFile(this->FileName, NcFile::ReadOnly);
    if(this->PointsFile->is_valid() == 0)
      {
      vtkErrorMacro(<< "Can't read file " << this->FileName);
      delete this->PointsFile;
      this->PointsFile = NULL;
      return 0;
      }
    this->SetCurrentFileName(this->FileName);
    }
  NcDim* timeDimension = this->PointsFile->get_dim("time");
  if(timeDimension == NULL)
    {
    vtkErrorMacro("Cannot find the number of time steps (time dimension).");
    return 0;
    }
  this->NumberOfTimeSteps = timeDimension->size();
217
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
218 219 220

  if (this->NumberOfTimeSteps > 0)
    {
221
    delete []this->TimeSteps;
222 223 224 225 226
    this->TimeSteps = new double[this->NumberOfTimeSteps];
    NcVar* timeVar = this->PointsFile->get_var("time");
    timeVar->get(this->TimeSteps, this->NumberOfTimeSteps);

    // Tell the pipeline what steps are available
227
    outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
228 229 230 231 232
              this->TimeSteps, this->NumberOfTimeSteps);

    // Range is required to get GUI to show things
    double tRange[2] = {this->TimeSteps[0],
                        this->TimeSteps[this->NumberOfTimeSteps - 1]};
233
    outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), tRange, 2);
234 235 236
    }
  else
    {
237 238
    outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
    outInfo->Remove(vtkStreamingDemandDrivenPipeline::TIME_RANGE());
239 240
    }

241
  outInfo->Set(
Berk Geveci's avatar
Berk Geveci committed
242
    CAN_HANDLE_PIECE_REQUEST(), 1);
243

244 245 246
  return 1;
}

247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
//----------------------------------------------------------------------------
int vtkNetCDFCAMReader::RequestUpdateExtent(
  vtkInformation *,
  vtkInformationVector **,
  vtkInformationVector *outputVector)
{
  if(this->FileName == NULL || this->ConnectivityFileName == NULL)
    {
    vtkWarningMacro("Missing a file name.");
    return 0;
    }
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  int piece =
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  int numPieces =
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());

  // make sure piece is valid
  if (piece < 0 || piece >= numPieces)
    {
268
    return 0;
269 270 271 272 273 274 275
    }

  return 1;
}

//----------------------------------------------------------------------------
int vtkNetCDFCAMReader::RequestData(
276
  vtkInformation *,vtkInformationVector **,vtkInformationVector *outputVector)
277 278 279 280 281 282 283 284 285 286 287 288
{
  if(this->FileName == NULL || this->ConnectivityFileName == NULL)
    {
    vtkWarningMacro("Missing a file name.");
    return 0;
    }

  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDebugMacro(<<"Reading NetCDF CAM file.");
289
  this->SetProgress(0);
290 291
  if(this->CurrentConnectivityFileName != NULL &&
     strcmp(this->CurrentConnectivityFileName, this->ConnectivityFileName) != 0)
292
    {
293 294 295
    delete this->ConnectivityFile;
    this->ConnectivityFile = NULL;
    this->SetCurrentConnectivityFileName(NULL);
296 297 298 299 300 301 302 303 304 305 306 307
    }
  if(this->ConnectivityFile == NULL)
    {
    this->ConnectivityFile = new NcFile(this->ConnectivityFileName,
                                        NcFile::ReadOnly);
    if(this->ConnectivityFile->is_valid() == 0)
      {
      vtkErrorMacro(<< "Can't read file " << this->ConnectivityFileName);
      delete this->ConnectivityFile;
      this->ConnectivityFile = NULL;
      return 0;
      }
308
    this->SetCurrentConnectivityFileName(this->ConnectivityFileName);
309 310
    }

311 312 313 314 315
  // Set the NetCDF error handler to not kill the application.
  // Upon exiting this method the error handler will be restored
  // to its previous state.
  NcError ncError(NcError::verbose_nonfatal);

316
  // read in the points first
317 318 319 320 321 322
  NcDim* levelsDimension = this->PointsFile->get_dim("lev");
  if(levelsDimension == NULL)
    {
    vtkErrorMacro("Cannot find the number of levels (lev dimension).");
    return 0;
    }
323
  long numCellLevels = levelsDimension->size() - 1;
324 325 326 327 328 329 330
  NcVar* levelsVar = this->PointsFile->get_var("lev");
  if(levelsVar == NULL)
    {
    vtkErrorMacro("Cannot find the number of levels (lev variable).");
    return 0;
    }
  if(levelsVar->num_dims() != 1 ||
331
     levelsVar->get_dim(0)->size() != numCellLevels+1)
332 333 334 335
    {
    vtkErrorMacro("The lev variable is not consistent.");
    return 0;
    }
336 337
  if(this->SingleLevel != 0)
    {
338
    numCellLevels = 1;
339
    }
340
  NcDim* dimension = this->PointsFile->get_dim("ncol");
341 342 343 344 345
  if(dimension == NULL)
    {
    vtkErrorMacro("Cannot find the number of points (ncol dimension).");
    return 0;
    }
346 347
  NcVar* lon = this->PointsFile->get_var("lon");
  NcVar* lat = this->PointsFile->get_var("lat");
348
  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
349 350
  output->SetPoints(points);

351
  long numFilePoints = dimension->size();
352 353
  if(lat == NULL || lon == NULL)
    {
354
    vtkErrorMacro("Cannot find coordinates (lat or lon variable).");
355 356 357 358 359
    return 0;
    }
  if(lat->type() == ncDouble)
    {
    points->SetDataTypeToDouble();
360 361 362
    points->SetNumberOfPoints(numFilePoints);
    std::vector<double> array(numFilePoints*2);
    if(!lon->get(&array[0], numFilePoints))
363 364 365
      {
      return 0;
      }
366
    if(!lat->get(&array[numFilePoints], numFilePoints))
367 368 369
      {
      return 0;
      }
370
    for(long i=0;i<numFilePoints;i++)
371
      {
372
      points->SetPoint(i, array[i], array[i+numFilePoints], numCellLevels+1);
373 374 375 376 377
      }
    }
  else
    {
    points->SetDataTypeToFloat();
378 379 380
    points->SetNumberOfPoints(numFilePoints);
    std::vector<float> array(numFilePoints*2);
    if(!lon->get(&array[0], numFilePoints))
381 382 383
      {
      return 0;
      }
384
    if(!lat->get(&array[numFilePoints], numFilePoints))
385 386 387
      {
      return 0;
      }
388
    for(long i=0;i<numFilePoints;i++)
389
      {
390
      points->SetPoint(i, array[i], array[i+numFilePoints], numCellLevels+1);
391 392
      }
    }
393 394
  this->SetProgress(.25);  // educated guess for progress

395 396 397 398
  // now read in the cell connectivity.  note that this is a periodic
  // domain and only the points on the left boundary are included in
  // the points file.  if a cell uses a point that is on the left
  // boundary and it should be on the right boundary we will have
399 400 401 402 403 404 405 406 407 408 409 410 411
  // to create that point.  That's what boundaryPoints is used for.
  // The (index + numFilePoints) gives us the new point id, and the value
  // for that in this array will correspond to the original point id that the
  // boundaryPoint is duplicate of.
  std::vector<vtkIdType> boundaryPoints;

  // To avoid creating multiple duplicates, we we a
  // vtkIncrementalOctreePointLocator.
  vtkSmartPointer<vtkIncrementalOctreePointLocator> locator =
    vtkSmartPointer<vtkIncrementalOctreePointLocator>::New();
  locator->SetDataSet(output); // dataset only has points right now.
  locator->BuildLocator();

412
  dimension = this->ConnectivityFile->get_dim("ncells");
413
  if(dimension == NULL)
414 415 416
    {
    vtkErrorMacro("Cannot find the number of cells (ncells dimension).");
    return 0;
417
    }
418 419
  NcVar* connectivity =
    this->ConnectivityFile->get_var("element_corners");
420 421
  if(connectivity == NULL)
    {
422
    vtkErrorMacro("Cannot find cell connectivity (element_corners dimension).");
423 424
    return 0;
    }
425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
  long numCellsPerLevel = dimension->size();


  int piece = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  int numPieces = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  int beginCellLevel, endCellLevel, beginCell, endCell;
  this->GetPartitioning(piece, numPieces, numCellLevels, numCellsPerLevel,
                        beginCellLevel, endCellLevel, beginCell, endCell);
  // the cells/levels assigned to this piece
  long numLocalCells = endCell-beginCell;
  int numLocalCellLevels = endCellLevel-beginCellLevel;
  std::vector<int> cellConnectivity(4*numLocalCells);
  connectivity->set_cur(0, beginCell);
  connectivity->get(&(cellConnectivity[0]), 4, numLocalCells);

  for(long i=0;i<numLocalCells;i++)
441 442
    {
    vtkIdType pointIds[4];
443
    double coords[4][3]; // assume quads here
444 445
    for(int j=0;j<4;j++)
      {
446
      pointIds[j] = cellConnectivity[i+j*numLocalCells]-1;
447
      points->GetPoint(pointIds[j], coords[j]);
448
      }
449
    if (IsCellInverted(coords) == true)
450
      {
451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488
      // First decide whether we're putting this cell on the 360 side (right) or on the
      // 0 side (left). We decide this based on which side will have the
      // smallest protrusion.
      double delta = 0.0;
      bool anchorLeft = false;
      for (int j=0; j < 4; ++j)
        {
        // We're assured that coords[j][0] is in the range [0, 360].
        // We just that fact to avoid having to do a std::abs() here.
        double rightDelta = (360.0 - coords[j][0]);
        double leftDelta = coords[j][0]; // i.e. (coords[j][0] - 0.0).
        if (IsZero(rightDelta) || IsZero(leftDelta) || rightDelta == leftDelta)
          {
          // if the point is equidistant from both ends or is one of the ends,
          // we let the other points in this cell dictate where the cell should
          // anchor since this point can easily be anchored on either side with
          // no side effects.
          continue;
          }
        if (rightDelta < leftDelta)
          {
          if (rightDelta > delta)
            {
            delta = rightDelta;
            anchorLeft = false;
            }
          }
        else
          {
          if (leftDelta > delta)
            {
            delta = leftDelta;
            anchorLeft = true;
            }
          }
        }
      // Once we've decided where we're anchoring we adjust the points.
      for (int j=0; j < 4; ++j)
489
        {
490
        if (anchorLeft)
491
          {
492 493 494 495
          // if coords[j] is closer to right (360), move it to the left.
          if ( (360.0 - coords[j][0]) < coords[j][0] )
            {
            coords[j][0] -= 360.0;
496 497
            }
          else
498 499
            {
            continue;
500 501
            }
          }
502
        else
503
          {
504 505 506 507
          // if coords[j] is closer to left (0), move it to the right
          if ( coords[j][0] < (360.0 - coords[j][0]) )
            {
            coords[j][0] += 360.0;
508 509
            }
          else
510 511
            {
            continue;
512 513
            }
          }
514 515 516 517 518 519 520 521 522 523 524 525
        // Okay, we have moved the coords. Update the boundaryPoints so which
        // original point id is this new point id a clone of.
        vtkIdType newPtId;
        if (locator->InsertUniquePoint(coords[j], newPtId) == 1)
          {
          // if a new point was indeed inserted, we need to update the
          // boundaryPoints to keep track of it.
          assert(newPtId >= numFilePoints && pointIds[j] < newPtId);
          assert(static_cast<vtkIdType>(boundaryPoints.size()) == (newPtId-numFilePoints));
          boundaryPoints.push_back(pointIds[j]);
          }
        cellConnectivity[i+j*numLocalCells] = (newPtId + 1); // note: 1-indexed.
526 527
        }
      }
528
    }
529
  locator = NULL; // release the locator memory.
530 531 532 533

  // we now have all of the points at a single level.  build them up
  // for the rest of the levels before creating the cells.
  vtkIdType numPointsPerLevel = points->GetNumberOfPoints();
534
  if(!this->SingleLevel)
535
    {
536
    // a hacky way to resize the points array without resetting the data
537
    points->InsertPoint(numPointsPerLevel*(numLocalCellLevels+1)-1, 0, 0, 0);
538
    for(vtkIdType pt=0;pt<numPointsPerLevel;pt++)
539
      {
540 541
      double point[3];
      points->GetPoint(pt, point);
542 543 544
      // need to start at 0 here since for multiple process the first
      // level will need to be replaced
      for(long lev=0;lev<numLocalCellLevels+1;lev++)
545
        {
546
        point[2] = numCellLevels - lev - beginCellLevel;
547 548
        points->SetPoint(pt+lev*numPointsPerLevel, point);
        }
549 550
      }
    }
551

552
  points->Modified();
553
  points->Squeeze();
554

555 556
  this->SetProgress(.5);  // educated guess for progress

557
  // Collect the time step requested
558 559 560
  vtkInformationDoubleKey* timeKey =
    static_cast<vtkInformationDoubleKey*>
    (vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
561 562 563 564

  double dTime = 0.0;
  if (outInfo->Has(timeKey))
    {
565
    dTime = outInfo->Get(timeKey);
566 567 568
    }

  // Actual time for the time step
569
  output->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(), dTime);
570 571 572 573 574 575 576 577 578

  // Index of the time step to request
  int timeStep = 0;
  while (timeStep < this->NumberOfTimeSteps &&
         this->TimeSteps[timeStep] < dTime)
    {
    timeStep++;
    }

579 580 581
  // now that we have the full set of points, read in any point
  // data with dimensions (time, lev, ncol) but read them in
  // by chunks of ncol since it will be a pretty big chunk of
582
  // memory that we'll have to break up anyways.
583 584 585
  for(int i=0;i<this->PointsFile->num_vars();i++)
    {
    NcVar* variable = this->PointsFile->get_var(i);
586 587 588 589
    if(this->SingleLevel == 0 && (variable->num_dims() != 3 ||
                                  strcmp(variable->get_dim(0)->name(), "time") != 0 ||
                                  strcmp(variable->get_dim(1)->name(), "lev") != 0 ||
                                  strcmp(variable->get_dim(2)->name(), "ncol") != 0) )
590
      { // not a 3D field variable
591 592
      continue;
      }
593 594 595
    else if(this->SingleLevel == 1 && ((variable->num_dims() != 2 ||
                                        strcmp(variable->get_dim(0)->name(), "time") != 0 ||
                                        strcmp(variable->get_dim(1)->name(), "ncol") != 0) ) )
596 597 598 599
      { // not a 2D field variable
      continue;
      }

600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
    vtkDoubleArray* doubleArray = NULL;
    vtkFloatArray* floatArray = NULL;
    if(variable->type() == ncDouble)
      {
      doubleArray = vtkDoubleArray::New();
      doubleArray->SetNumberOfTuples(points->GetNumberOfPoints());
      doubleArray->SetName(variable->name());
      output->GetPointData()->AddArray(doubleArray);
      doubleArray->Delete();
      }
    else
      {
      floatArray = vtkFloatArray::New();
      floatArray->SetNumberOfTuples(points->GetNumberOfPoints());
      floatArray->SetName(variable->name());
      output->GetPointData()->AddArray(floatArray);
      floatArray->Delete();
      }
618
    if(this->SingleLevel == 0)
619
      {
620
      for(long lev=0;lev<numLocalCellLevels+1;lev++)
621
        {
622
        variable->set_cur(timeStep, lev+beginCellLevel, 0);
623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642
        if(doubleArray)
          {
          if(!variable->get(doubleArray->GetPointer(0)+lev*numPointsPerLevel,
                            1, 1, numFilePoints))
            {
            vtkErrorMacro("Problem getting NetCDF variable " << variable->name());
            return 0;
            }
          }
        else
          {
          if(!variable->get(floatArray->GetPointer(0)+lev*numPointsPerLevel,
                            1, 1, numFilePoints))
            {
            vtkErrorMacro("Problem getting NetCDF variable " << variable->name());
            return 0;
            }
          }
        }
      }
643
    else if(this->SingleLevel == 1)
644
      {
645
      variable->set_cur(timeStep, 0);
646 647
      if(doubleArray)
        {
648
        if(!variable->get(doubleArray->GetPointer(0), 1, numFilePoints))
649
          {
650
          vtkErrorMacro("Problem getting NetCDF variable " << variable->name());
651 652 653 654 655
          return 0;
          }
        }
      else
        {
656
        if(!variable->get(floatArray->GetPointer(0), 1, numFilePoints))
657
          {
658
          vtkErrorMacro("Problem getting NetCDF variable " << variable->name());
659 660 661 662
          return 0;
          }
        }
      }
663
    }
664 665

  // we have to copy the values from the left side to the right side
666
  vtkPointData* pointData = output->GetPointData();
667 668 669
  pointData->CopyAllOn();
  pointData->CopyAllocate(output->GetPointData(),
                          output->GetNumberOfPoints());
670

671
  vtkIdType newPtId=0;
672 673
  for (std::vector<vtkIdType>::const_iterator it=
        boundaryPoints.begin(); it!=boundaryPoints.end(); ++it, ++newPtId)
674
    {
675
    for(long lev=0;lev<numLocalCellLevels+1-this->SingleLevel;lev++)
676
      {
677 678 679
      vtkIdType srcId = (*it) + lev * numPointsPerLevel;
      vtkIdType destId = (newPtId + numFilePoints) + lev * numPointsPerLevel;
      pointData->CopyData(pointData, srcId, destId);
680 681
      }
    }
682

683
  // add in level data for each plane which corresponds to an average pressure
684 685 686 687 688 689 690 691 692 693
  // if we are loading a volumetric grid
  if(this->SingleLevel == 0)
    {
    std::vector<float> levelData(numLocalCellLevels+1);
    levelsVar->set_cur(beginCellLevel);
    levelsVar->get(&levelData[0], numLocalCellLevels+1);
    vtkNew<vtkFloatArray> levelPointData;
    levelPointData->SetName(levelsVar->name());
    levelPointData->SetNumberOfTuples(points->GetNumberOfPoints());
    for(long j=0;j<numLocalCellLevels+1;j++)
694
      {
695 696 697 698
      for(vtkIdType i=0;i<numPointsPerLevel;i++)
        {
        levelPointData->SetValue(j*numPointsPerLevel+i, levelData[j]);
        }
699
      }
700
    output->GetPointData()->AddArray(levelPointData.GetPointer());
701
    }
702

703 704
  this->SetProgress(.75);  // educated guess for progress

705
  // now we actually create the cells
706 707 708 709 710 711 712 713 714
  if(this->SingleLevel == 1)
    {
    output->Allocate(numLocalCells);
    }
  else
    {
    output->Allocate(numLocalCells*numLocalCellLevels);
    }
  for(long i=0;i<numLocalCells;i++)
715 716 717 718
    {
    vtkIdType pointIds[4];
    for(int j=0;j<4;j++)
      {
719
      pointIds[j] = cellConnectivity[i+j*numLocalCells]-1;
720
      }
721
    if(this->SingleLevel == 0)
722
      { // volumetric grid
723
      for(int lev=0;lev<numLocalCellLevels;lev++)
724
        {
725 726
        vtkIdType hexIds[8];
        for(int j=0;j<4;j++)
727
          {
728 729
          hexIds[j] = pointIds[j]+lev*numPointsPerLevel;
          hexIds[j+4] = pointIds[j]+(1+lev)*numPointsPerLevel;
730
          }
731
        output->InsertNextCell(VTK_HEXAHEDRON, 8, hexIds);
732 733
        }
      }
734
    else if(this->SingleLevel == 1)
735 736
      { // surface grid
      output->InsertNextCell(VTK_QUAD, 4, pointIds);
737
      }
738 739
    }

740 741 742 743 744 745 746
  if(numLocalCells != numCellsPerLevel)
    {
    // we have extra points that are not connected to any cells
    //vtkNew<vtkCleanUnstructuredGrid> cleanGrid;
    //cleanGrid->SetInput(output);
    }

747 748
  vtkDebugMacro(<<"Read " << output->GetNumberOfPoints() <<" points,"
                << output->GetNumberOfCells() <<" cells.\n");
749 750 751 752

  return 1;
}

753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
//----------------------------------------------------------------------------
bool vtkNetCDFCAMReader::GetPartitioning(
  int piece, int numPieces,int numLevels, int numCellsPerLevel,
  int & beginLevel, int & endLevel, int & beginCell, int & endCell)
{
  // probably not the best way to partition the data but should
  // be sufficient for development.
  if(numPieces <= 0 || piece < 0 || piece >= numPieces)
    {
    vtkErrorMacro("Bad piece information for partitioning.");
    return false;
    }
  if(numPieces == 1)
    {
    beginLevel = 0;
    endLevel = numLevels;
    beginCell = 0;
    endCell = numCellsPerLevel;
    return true;
    }
  if(numPieces <= numLevels)
    {
    beginLevel = piece*numLevels/numPieces;
    endLevel = (piece+1)*numLevels/numPieces;
    beginCell = 0;
    endCell = numCellsPerLevel;
    return true;
    }

  int levelsPerPiece = vtkMath::Ceil(numLevels/static_cast<double>(numPieces));
  int piecesPerLevel = vtkMath::Ceil(numPieces/static_cast<double>(numLevels));
  int numOverworkedPieces = piecesPerLevel/levelsPerPiece*numLevels - numPieces;
  bool evenOverworked = (piecesPerLevel % 2 == 0 || numOverworkedPieces == 0);
  if(piece < numOverworkedPieces)
    {
    if(evenOverworked)
      {
      beginLevel = 2*piece/piecesPerLevel;
      endLevel = beginLevel + 1;
      int remainder = piece % (piecesPerLevel/2);
      beginCell = remainder * numCellsPerLevel * 2 / piecesPerLevel;
      endCell = (remainder + 1)* numCellsPerLevel * 2 / piecesPerLevel;
      }
    else
      {
      beginLevel = 2*piece/(piecesPerLevel-1);
      endLevel = beginLevel + 1;
      int remainder = piece % ((piecesPerLevel-1)/2);
      beginCell = remainder * numCellsPerLevel * 2 / piecesPerLevel;
      endCell = (remainder + 1)* numCellsPerLevel * 2 / piecesPerLevel;
      }
    }
  else // underworked pieces
    {
    if( evenOverworked == false && piece - numOverworkedPieces < 2*numOverworkedPieces/(piecesPerLevel-1) )
      { // fillers for levels that also have overworked pieces working on them
      beginLevel = piece - numOverworkedPieces;
      beginCell = numCellsPerLevel*(piecesPerLevel-1)/piecesPerLevel;
      endCell = numCellsPerLevel;
      }
    else
      {
      int fakePiece = numOverworkedPieces+piece; // take into account overworked pieces
      beginLevel = fakePiece / piecesPerLevel;
      int remainder = fakePiece % piecesPerLevel;
      beginCell = remainder * numCellsPerLevel / piecesPerLevel;
      endCell = (remainder + 1)*numCellsPerLevel / piecesPerLevel;
      }
    endLevel = beginLevel + 1;
    }

  return true;
}

827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843
//----------------------------------------------------------------------------
#if !defined(VTK_LEGACY_REMOVE)
void vtkNetCDFCAMReader::SetCellLayerRight(int)
{
  VTK_LEGACY_BODY(vtkNetCDFCAMReader::SetCellLayerRight, "VTK 6.3");
}
#endif

//----------------------------------------------------------------------------
#if !defined(VTK_LEGACY_REMOVE)
int vtkNetCDFCAMReader::GetCellLayerRight()
{
  VTK_LEGACY_BODY(vtkNetCDFCAMReader::GetCellLayerRight, "VTK 6.3");
  return 0;
}
#endif

844 845 846 847
//----------------------------------------------------------------------------
void vtkNetCDFCAMReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
848 849 850 851 852
  os << indent << "FileName: "
     << (this->FileName ? this->FileName : "(NULL)") << endl;
  os << indent << "ConnectivityFileName: " <<
    (this->ConnectivityFileName ? this->ConnectivityFileName : "(NULL)")
     << endl;
853
  os << indent << "SingleLevel: " << this->SingleLevel << endl;
854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869
  if(this->PointsFile)
    {
    os << indent << "PointsFile: " << this->PointsFile << endl;
    }
  else
    {
    os << indent << "PointsFile: (NULL)" << endl;
    }
  if(this->ConnectivityFile)
    {
    os << indent << "ConnectivityFile: " << this->ConnectivityFile << endl;
    }
  else
    {
    os << indent << "ConnectivityFile: (NULL)" << endl;
    }
870
}