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
}