vtkXdmfReader.cxx 91.1 KB
Newer Older
1
/*=========================================================================
Jerry Clarke's avatar
Jerry Clarke committed
2

3
4
5
6
7
  Program:   Visualization Toolkit
  Module:    vtkXdmfReader.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$
Jerry Clarke's avatar
Jerry Clarke committed
8
9


Berk Geveci's avatar
Berk Geveci committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
  Copyright (c) 1993-2001 Ken Martin, Will Schroeder, Bill Lorensen  
  All rights reserved.

  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that the following conditions are met:

  * Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.

  * Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

  * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
  of any contributors may be used to endorse or promote products derived
  from this software without specific prior written permission.

  * Modified source versions must be plainly marked as such, and must not be
  misrepresented as being the original software.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

  =========================================================================*/
Andy Cedilnik's avatar
Andy Cedilnik committed
42
#include "vtkXdmfReader.h"
43
44
45
46
47

#include "vtkCallbackCommand.h"
#include "vtkDataArraySelection.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
48
49
#include "vtkInformation.h"
#include "vtkInformationVector.h"
50
51
#include "vtkObjectFactory.h"
#include "vtkRectilinearGrid.h"
52
#include "vtkStreamingDemandDrivenPipeline.h"
53
54
55
#include "vtkStructuredGrid.h"
#include "vtkUnstructuredGrid.h"
#include "vtkXdmfDataArray.h"
56
57
58
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkCellArray.h"
59
#include "vtkCharArray.h"
60
#include "vtkXMLParser.h"
61
#include "vtkImageData.h"
62
#include "vtkUniformGrid.h"
63
#include "vtkMultiGroupDataInformation.h"
64
#include "vtkMultiGroupDataSet.h"
65
66
#include "vtkCompositeDataPipeline.h"
#include "vtkMultiProcessController.h"
67

68
#include "XdmfArray.h"
Jerry Clarke's avatar
Jerry Clarke committed
69
#include "XdmfAttribute.h"
70
71
#include "XdmfDOM.h"
#include "XdmfDataDesc.h"
Jerry Clarke's avatar
Jerry Clarke committed
72
#include "XdmfDataItem.h"
73
#include "XdmfGrid.h"
Jerry Clarke's avatar
Jerry Clarke committed
74
#include "XdmfTopology.h"
Andy Cedilnik's avatar
Andy Cedilnik committed
75
#include "XdmfGeometry.h"
Jerry Clarke's avatar
Jerry Clarke committed
76

77
#include <sys/stat.h>
78
79
#include <vtkstd/set>
#include <vtkstd/map>
80
81
#include <vtkstd/string>
#include <vtkstd/vector>
82
83
84
#include <assert.h>

#define USE_IMAGE_DATA // otherwise uniformgrid
Jerry Clarke's avatar
Jerry Clarke committed
85

Andy Cedilnik's avatar
Andy Cedilnik committed
86
vtkStandardNewMacro(vtkXdmfReader);
87
vtkCxxRevisionMacro(vtkXdmfReader, "1.27");
88
89

vtkCxxSetObjectMacro(vtkXdmfReader,Controller,vtkMultiProcessController);
90
91
92
93
94
95
96
97

#if defined(_WIN32) && (defined(_MSC_VER) || defined(__BORLANDC__))
#  include <direct.h>
#  define GETCWD _getcwd
#else
#include <unistd.h>
#  define GETCWD getcwd
#endif
Jerry Clarke's avatar
Jerry Clarke committed
98

99
100
101
#define vtkMAX(x, y) (((x)>(y))?(x):(y))
#define vtkMIN(x, y) (((x)<(y))?(x):(y))

102
103
#define PRINT_EXTENT(x) "[" << (x)[0] << " " << (x)[1] << " " << (x)[2] << " " << (x)[3] << " " << (x)[4] << " " << (x)[5] << "]" 

104
//============================================================================
105
106
107
108
109
110
class vtkXdmfReaderGrid
{
public:
  vtkXdmfReaderGrid() : XMGrid(0), DataDescription(0) {}
  ~vtkXdmfReaderGrid()
    {
Berk Geveci's avatar
Berk Geveci committed
111
      delete this->XMGrid;
112
113
114
115
116
    }

  XdmfGrid       *XMGrid;
  XdmfDataDesc   *DataDescription;
  vtkstd::string Name;
117
  int Level;
118
119
};

120
//============================================================================
121
122
123
class vtkXdmfReaderGridCollection
{
public:
124
125
  vtkXdmfReaderGrid* GetXdmfGrid(const char* gridName,
                                 int level);
126
127
128

  typedef vtkstd::map<vtkstd::string, vtkXdmfReaderGrid*> SetOfGrids;
  SetOfGrids Grids;
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
  // Update number of levels and number of datasets.
  void UpdateCounts();
  
  // Return the last number of levels computed by UpdateCounts.
  int GetNumberOfLevels()
    {
      return this->NumberOfLevels;
    }
  // Return the last number of dataset for a given level
  // computed by UpdateCounts.
  int GetNumberOfDataSets(int level)
    {
      assert("pre: valid_level" && level>=0 && level<GetNumberOfLevels());
      return this->NumberOfDataSets[level];
    }
  
protected:
  int NumberOfLevels;
  vtkstd::vector<int> NumberOfDataSets;
148
149
150
};

//----------------------------------------------------------------------------
151
vtkXdmfReaderGrid* vtkXdmfReaderGridCollection::GetXdmfGrid(
152
153
  const char *gridName,
  int level)
154
{
155
  //find the grid with that name or make one up
156
157
158
159
160
  vtkXdmfReaderGridCollection::SetOfGrids::iterator it = this->Grids.find(gridName);
  if ( it == this->Grids.end() || it->second == 0 )
    {
    this->Grids[gridName] = new vtkXdmfReaderGrid;
    }
161
  //sets its level
162
  this->Grids[gridName]->Level=level;
163
164
165
  return this->Grids[gridName];
}

166
167
168
169
170
171
172
//----------------------------------------------------------------------------
void vtkXdmfReaderGridCollection::UpdateCounts()
{
  // Update the number of levels.
  vtkXdmfReaderGridCollection::SetOfGrids::iterator it;
  vtkXdmfReaderGridCollection::SetOfGrids::iterator itEnd;
  
173
  //find maximum level of any of my grids
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  int maxLevel=0;
  it=this->Grids.begin();
  itEnd=this->Grids.end();
  while(it!=itEnd)
    {
    int level=it->second->Level;
    if(level>maxLevel)
      {
      maxLevel=level;
      }
    ++it;
    }
  this->NumberOfLevels=maxLevel+1;
  this->NumberOfDataSets.resize(this->NumberOfLevels);
188
189

  //clear the number of datasets at each level array
190
191
192
193
194
195
196
197
  int l=this->NumberOfLevels;
  int i=0;
  while(i<l)
    {
    this->NumberOfDataSets[i]=0;
    ++i;
    }
  
198
  //count up the number of datasets at each level and save in the array
199
200
201
202
203
204
205
206
207
  it=this->Grids.begin();
  while(it!=itEnd)
    {
    int level=it->second->Level;
    ++this->NumberOfDataSets[level];
    ++it;
    }
}

208
//============================================================================
209
210
211
212
class vtkXdmfReaderActualGrid
{
public:
  vtkXdmfReaderActualGrid() : Enabled(0), Grid(0), Collection(0) {}
Dave Demarle's avatar
Dave Demarle committed
213
214
  ~vtkXdmfReaderActualGrid() 
  {
215
  /*
Dave Demarle's avatar
Dave Demarle committed
216
217
218
219
220
221
222
223
224
225
226
227
228
  if (this->Collection)
    {
    vtkXdmfReaderGridCollection::SetOfGrids::iterator it = this->Collection->Grids.begin();
    while ( it != this->Collection->Grids.end() )
      {
      if (it->second)
        {
        delete it->second;
        it->second = NULL;
        }
      it++;
      }
    };
229
  */
Dave Demarle's avatar
Dave Demarle committed
230
  }
231

232
233
234
235
236
  int Enabled;
  vtkXdmfReaderGrid* Grid;
  vtkXdmfReaderGridCollection* Collection;
};

237
//============================================================================
Andy Cedilnik's avatar
Andy Cedilnik committed
238
class vtkXdmfReaderInternal
239
240
{
public:
241
242
  vtkXdmfReaderInternal()
    {
Berk Geveci's avatar
Berk Geveci committed
243
244
245
      this->DataItem = 0;
      this->ArrayConverter = vtkXdmfDataArray::New();
      this->NumberOfEnabledActualGrids=0;
246
247
248
    }
  ~vtkXdmfReaderInternal()
    {
Berk Geveci's avatar
Berk Geveci committed
249
250
251
252
253
254
255
256
257
      if ( this->DataItem )
        {
        // cout << "....Deleting DataItem " << this->DataItem << endl;
        delete this->DataItem;
        // cout << "....Back from Deleting DataItem " << this->DataItem << endl;
        this->DataItem = 0;
        }
      this->ArrayConverter->Delete();
      this->ArrayConverter = 0;
258
259
    }

260
  typedef vtkstd::vector<vtkstd::string> StringListType;
261
262
  StringListType DomainList;
  XdmfXmlNode DomainPtr;
263
264

  int RequestSingleGridData(const char* currentGridName,
265
266
                            vtkXdmfReaderGrid *grid,
                            vtkInformation *outInfo,
267
268
                            vtkDataObject *output,
                            int isSubBlock);
269
270
271
272
  
  int RequestActualGridData(const char* currentGridName,
                            vtkXdmfReaderActualGrid* currentActualGrid,
                            int outputGrid,
273
                            int numberOfGrids,
274
275
276
277
278
279
280
281
                            vtkInformationVector *outputVector);
  
  // outInfo is null in the  multi-block case.
  int RequestSingleGridInformation(vtkXdmfReaderGrid *grid,
                                   vtkInformation *outInfo);

  int RequestActualGridInformation(vtkXdmfReaderActualGrid* currentActualGrid,
                                   int outputGrid,
Berk Geveci's avatar
Berk Geveci committed
282
                                   int numberOfGrids,
283
284
285
                                   vtkInformationVector* outputVector);
  
  typedef vtkstd::map<vtkstd::string,vtkXdmfReaderActualGrid> MapOfActualGrids;
286
  MapOfActualGrids ActualGrids;
287
288
  int NumberOfEnabledActualGrids;
  
289
290
291
  vtkXdmfReaderGridCollection* GetCollection(const char* collectionName);
  vtkXdmfReaderActualGrid* GetGrid(const char* gridName);
  vtkXdmfReaderActualGrid* GetGrid(int idx);
292
  vtkXdmfReaderGrid *GetXdmfGrid(const char *gridName,
293
294
                                 const char *collectionName,
                                 const char *levelName);
295
296

  vtkXdmfReader* Reader;
Jerry Clarke's avatar
Jerry Clarke committed
297
  XdmfDataItem *DataItem;
298
299
300

  // For converting arrays from XDMF to VTK format
  vtkXdmfDataArray *ArrayConverter;
301
302
};

303
//----------------------------------------------------------------------------
304
305
vtkXdmfReaderGridCollection* vtkXdmfReaderInternal::GetCollection(
  const char* collectionName)
306
307
308
309
310
{
  if ( !collectionName )
    {
    return 0;
    }
311

312
  vtkXdmfReaderActualGrid* actualGrid = &this->ActualGrids[collectionName];
313

314
315
316
317
318
319
320
321
322
  if ( !actualGrid->Collection )
    {
    if ( actualGrid->Grid )
      {
      cerr << "Trying to create collection with the same name as an existing grid" << endl;
      return 0;
      }
    actualGrid->Collection = new vtkXdmfReaderGridCollection;
    }
323

324
325
326
327
  return actualGrid->Collection;
}

//----------------------------------------------------------------------------
328
vtkXdmfReaderGrid* vtkXdmfReaderInternal::GetXdmfGrid(
329
330
331
  const char *gridName,
  const char *collectionName,
  const char *levelName)
332
333
334
335
336
337
338
339
{
  if ( !gridName )
    {
    return 0;
    }

  if ( collectionName )
    {
340
    vtkXdmfReaderGridCollection *collection=this->GetCollection(collectionName);
341
342
343
344
    if ( !collection )
      {
      return 0;
      }
345
346
347
348
349
350
351
    int level;
    if(levelName==0)
      {
      level=0;
      }
    else
      {
352
353
354
      char *tmp=new char[strlen(levelName)+1];
      strcpy(tmp,levelName);
      istrstream s(tmp,strlen(tmp));
355
356
357
      s>>level;
      if(level<0)
        {
358
359
        cerr << "Expect a positive Level value, level=" << level <<endl;
        delete[] tmp;
360
361
        return 0;
        }
362
363
364
365
      else
        {
        delete[] tmp;
        }
366
367
      }
    return collection->GetXdmfGrid(gridName,level);
368
369
370
371
372
373
374
375
    }

  vtkXdmfReaderActualGrid* grid = this->GetGrid(gridName);
  if ( grid->Collection )
    {
    cerr << "Trying to create a grid with the same name as an existing collection" << endl;
    return 0;
    }
376
  grid->Grid = new vtkXdmfReaderGrid; 
377
378
379
380
381
382
383
384
385
386
387
388
  return grid->Grid;
}

//----------------------------------------------------------------------------
vtkXdmfReaderActualGrid* vtkXdmfReaderInternal::GetGrid(const char* gridName)
{
  return &this->ActualGrids[gridName];
}

//----------------------------------------------------------------------------
vtkXdmfReaderActualGrid* vtkXdmfReaderInternal::GetGrid(int idx)
{
389
390
391
392
  if ( idx < 0 )
    {
    return 0;
    }
393
394
395
  vtkXdmfReaderInternal::MapOfActualGrids::iterator it;
  int cnt = 0;
  for ( it = this->ActualGrids.begin();
Berk Geveci's avatar
Berk Geveci committed
396
397
        it != this->ActualGrids.end();
        ++ it )
398
399
400
401
402
403
404
405
406
407
    {
    if ( cnt == idx )
      {
      return &it->second;
      }
    cnt ++;
    }
  return 0;
}

408
//----------------------------------------------------------------------------
409
410
411
412
413
414
int vtkXdmfReaderInternal::RequestActualGridData(
  const char* vtkNotUsed(currentGridName),
  vtkXdmfReaderActualGrid* currentActualGrid,
  int outputGrid,
  int vtkNotUsed(numberOfGrids),
  vtkInformationVector *outputVector)
415
{
416
417
  vtkInformation *info=outputVector->GetInformationObject(0);
  int procId=info->Get(
418
    vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()); 
419
420
  int nbProcs=info->Get(
    vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
421
  
422
423
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkMultiGroupDataSet *mgd=vtkMultiGroupDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
424

425
  if (!currentActualGrid->Collection)
426
    {
427
    return 0; //work around an undo/redo bug
Jerry Clarke's avatar
Jerry Clarke committed
428
    }
429

430
431
432
433
434
435
  unsigned int numberOfDataSets=currentActualGrid->Collection->Grids.size();
    
  currentActualGrid->Collection->UpdateCounts();
  int levels=currentActualGrid->Collection->GetNumberOfLevels();
  // mgd->SetNumberOfLevels(levels);
  // int levels = 1;
436
  
437
438
439
440
441
442
443
  int level=0;
  mgd->SetNumberOfDataSets(outputGrid, currentActualGrid->Collection->GetNumberOfDataSets(level));
  while(level<levels)
    {
    // mgd->SetNumberOfDataSets(level,currentActualGrid->Collection->GetNumberOfDataSets(level));
    ++level;
    }
444
  
445
446
447
448
449
450
451
452
453
454
  vtkXdmfReaderGridCollection::SetOfGrids::iterator gridIt;
  vtkXdmfReaderGridCollection::SetOfGrids::iterator gridItEnd;
    
  int numBlocksPerProcess=numberOfDataSets/nbProcs;
  int leftOverBlocks=numberOfDataSets-(numBlocksPerProcess*nbProcs);
    
  int blockStart;
  int blockEnd;
    
  if(procId<leftOverBlocks)
Andy Cedilnik's avatar
Andy Cedilnik committed
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
489
    blockStart=(numBlocksPerProcess+1)*procId;
    blockEnd=blockStart+(numBlocksPerProcess+1)-1;
    }
  else
    {
    blockStart=numBlocksPerProcess*procId+leftOverBlocks;
    blockEnd=blockStart+numBlocksPerProcess-1;
    }
    
  gridIt=currentActualGrid->Collection->Grids.begin();
  gridItEnd=currentActualGrid->Collection->Grids.end();
  int result=1;
    
  int datasetIdx=0;
    
  vtkMultiGroupDataInformation *compInfo =
    vtkMultiGroupDataInformation::SafeDownCast(
      info->Get(vtkCompositeDataPipeline::COMPOSITE_DATA_INFORMATION()));
    
  // currentIndex in each level.
  vtkstd::vector<int> currentIndex(levels);
  level=0;
  while(level<levels)
    {
    currentIndex[level]=0;
    ++level;
    }
    
  while(gridIt != gridItEnd && result)
    {
    vtkXdmfReaderGrid *subgrid=gridIt->second;
    level=subgrid->Level;
    int index=currentIndex[level];
    if(datasetIdx<blockStart || datasetIdx>blockEnd)
490
      {
491
492
      // mgd->SetDataSet(level,index,0); // empty, on another processor
      mgd->SetDataSet(outputGrid, index, 0); // empty, on another processor
493
      }
494
    else
495
      {
496
497
      XdmfGrid *xdmfGrid=subgrid->XMGrid;
      if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
498
        {
499
500
501
502
503
504
505
506
507
508
509
510
511
        vtkUnstructuredGrid *ds=vtkUnstructuredGrid::New();
        ds->SetMaximumNumberOfPieces(1);
        // mgd->SetDataSet(level,index,ds);
        mgd->SetDataSet(outputGrid, index, ds);
        ds->Delete();
        } 
      else if( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DSMESH ||
               xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DSMESH )
        {
        vtkStructuredGrid *ds=vtkStructuredGrid::New();
        // mgd->SetDataSet(level,index,ds);
        mgd->SetDataSet(outputGrid, index, ds);
        ds->Delete();
512
        }
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
      else if 
        (xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DCORECTMESH ||
         xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DCORECTMESH )
        {
#ifdef USE_IMAGE_DATA
        vtkImageData *ds=vtkImageData::New();
#else
        vtkUniformGrid *ds=vtkUniformGrid::New();
#endif
        //mgd->SetDataSet(level,index,ds);
        mgd->SetDataSet(outputGrid,index,ds);
        ds->Delete();
        }
      else if ( 
        xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DRECTMESH ||
        xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DRECTMESH )
        {
        vtkRectilinearGrid *ds=vtkRectilinearGrid::New();
        //mgd->SetDataSet(level,index,ds);
        mgd->SetDataSet(outputGrid,index,ds);
        ds->Delete();
        }
      else
        {
        // Unknown type for this sub grid. 
        return 0;
        }
      vtkDataObject *ds=mgd->GetDataSet(outputGrid,index);
541
      vtkInformation *subInfo=compInfo->GetInformation(outputGrid,index); 
542
      result=this->RequestSingleGridData("",gridIt->second,subInfo,ds,1);
543
      }
544
545
546
547
    ++currentIndex[level];
    ++gridIt;
    ++datasetIdx;
    this->Reader->UpdateProgress(1.0 * datasetIdx / numberOfDataSets);
Andy Cedilnik's avatar
Andy Cedilnik committed
548
    }
549
  return result;
Jerry Clarke's avatar
Jerry Clarke committed
550
551
}

552
//----------------------------------------------------------------------------
553
554
555
556
557
558
int vtkXdmfReaderInternal::RequestSingleGridData(
  const char* currentGridName,
  vtkXdmfReaderGrid *grid,
  vtkInformation *outInfo,
  vtkDataObject *output,
  int isSubBlock)
559
{
560
561
562
563
564
565
566
567
568
569
570
571
  int *readerStride = this->Reader->GetStride();
  
  vtkDataArraySelection* pointDataArraySelection = 
    this->Reader->GetPointDataArraySelection();
  vtkDataArraySelection* cellDataArraySelection = 
    this->Reader->GetCellDataArraySelection();
  
  // Handle single grid
  XdmfGrid* xdmfGrid = grid->XMGrid;
  XdmfDOM* xdmfDOM = xdmfGrid->GetDOM();
  
  if ( !grid->DataDescription )
572
    {
573
574
    grid->DataDescription = xdmfGrid->GetTopology()->GetShapeDesc();
    //continue;
575
576
    }
  
577
578
579
580
581
582
583
584
  vtkDebugWithObjectMacro(this->Reader, 
                          << "Reading Heavy Data for " 
                          << xdmfGrid->GetName());
  xdmfGrid->Update();
  
  // True for all 3d datasets except unstructured grids
  int globalrank = 3;
  switch(xdmfGrid->GetTopology()->GetTopologyType())
585
    {
586
587
    case XDMF_2DSMESH: case XDMF_2DRECTMESH: case XDMF_2DCORECTMESH:
      globalrank = 2;
588
    }
Jerry Clarke's avatar
Jerry Clarke committed
589
  if ( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED )
590
591
592
593
594
595
    {
    globalrank = 1;
    }
  
  
  vtkIdType cc;
Jerry Clarke's avatar
Jerry Clarke committed
596
597
  XdmfXmlNode attrNode;
  XdmfXmlNode dataNode;
598
599
600
601
602
603
604
605
  XdmfInt64 start[4]  = { 0, 0, 0, 0 };
  XdmfInt64 stride[4] = { 1, 1, 1, 1 };
  XdmfInt64 count[4] = { 0, 0, 0, 0 };
  grid->DataDescription->GetShape(count);

  int upext[6];
  int whext[6];
  
Jerry Clarke's avatar
Jerry Clarke committed
606
  if( xdmfGrid->GetTopology()->GetClass() != XDMF_UNSTRUCTURED)
607
    {
608
    if(isSubBlock)
609
      {
610
611
612
      // the composite pipeline does not set update extent, so just take
      // the whole extent for as update extent.
      outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), upext);
613
      }
614
    else
615
      {
616
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), upext);
617
      }
618
    
619
620
    outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), whext);
    
621
622
623
    start[2] = vtkMAX(0, upext[0]);
    start[1] = vtkMAX(0, upext[2]);
    start[0] = vtkMAX(0, upext[4]);
624
    
625
626
627
    count[2] = upext[1] - upext[0];
    count[1] = upext[3] - upext[2];
    count[0] = upext[5] - upext[4];
628
    }
629
  
630
  XdmfGeometry  *Geometry = xdmfGrid->GetGeometry();
631
632
  
  
633
  // Read Topology for Unstructured Grid
Jerry Clarke's avatar
Jerry Clarke committed
634
  if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
635
636
637
638
639
640
641
642
643
644
645
    {
    vtkUnstructuredGrid  *vGrid = static_cast<vtkUnstructuredGrid *>(output);
    vtkCellArray                *verts;
    XdmfInt32           vType;
    XdmfInt32           NodesPerElement;
    vtkIdType           NumberOfElements;
    vtkIdType           i, j, index;    
    XdmfInt64           Length, *Connections;
    vtkIdType           *connections;
    int                 *cell_types, *ctp;
    
Berk Geveci's avatar
Berk Geveci committed
646
647
648
    vtkDebugWithObjectMacro(
      this->Reader,<< "Unstructured Topology is " 
      << xdmfGrid->GetTopology()->GetTopologyTypeAsString());
Jerry Clarke's avatar
Jerry Clarke committed
649
    switch ( xdmfGrid->GetTopology()->GetTopologyType() )
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
      {
      case  XDMF_POLYVERTEX :
        vType = VTK_POLY_VERTEX;
        break;
      case  XDMF_POLYLINE :
        vType = VTK_POLY_LINE;
        break;
      case  XDMF_POLYGON :
        vType = VTK_POLYGON;
        break;
      case  XDMF_TRI :
        vType = VTK_TRIANGLE;
        break;
      case  XDMF_QUAD :
        vType = VTK_QUAD;
        break;
      case  XDMF_TET :
        vType = VTK_TETRA;
        break;
      case  XDMF_PYRAMID :
        vType = VTK_PYRAMID;
Jerry Clarke's avatar
Jerry Clarke committed
671
        break;
672
673
674
675
676
677
      case  XDMF_WEDGE :
        vType = VTK_WEDGE;
        break;
      case  XDMF_HEX :
        vType = VTK_HEXAHEDRON;
        break;
Jerry Clarke's avatar
Jerry Clarke committed
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
      case  XDMF_EDGE_3 :
        vType = VTK_QUADRATIC_EDGE ;
        break;
      case  XDMF_TRI_6 :
        vType = VTK_QUADRATIC_TRIANGLE ;
        break;
      case  XDMF_QUAD_8 :
        vType = VTK_QUADRATIC_QUAD ;
        break;
      case  XDMF_TET_10 :
        vType = VTK_QUADRATIC_TETRA ;
        break;
      case  XDMF_PYRAMID_13 :
        vType = VTK_QUADRATIC_PYRAMID ;
        break;
      case  XDMF_WEDGE_15 :
        vType = VTK_QUADRATIC_WEDGE ;
        break;
      case  XDMF_HEX_20 :
        vType = VTK_QUADRATIC_HEXAHEDRON ;
        break;
Jerry Clarke's avatar
Jerry Clarke committed
699
700
701
      case XDMF_MIXED :
        vType = -1;
        break;
702
      default :
Berk Geveci's avatar
Berk Geveci committed
703
704
        XdmfErrorMessage("Unknown Topology Type = " 
                         << xdmfGrid->GetTopology()->GetTopologyType());
705
        return 1;
706
      }
707
    if( xdmfGrid->GetTopology()->GetTopologyType() != XDMF_MIXED)
Berk Geveci's avatar
Berk Geveci committed
708
      {
709
710
      NodesPerElement = xdmfGrid->GetTopology()->GetNodesPerElement();
      if ( xdmfGrid->GetTopology()->GetConnectivity()->GetRank() == 2 )
Berk Geveci's avatar
Berk Geveci committed
711
        {
712
713
        NodesPerElement = 
          xdmfGrid->GetTopology()->GetConnectivity()->GetDimension(1);
Berk Geveci's avatar
Berk Geveci committed
714
        }
Jerry Clarke's avatar
Jerry Clarke committed
715
    
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
      /* Create Cell Type Array */
      Length = xdmfGrid->GetTopology()->GetConnectivity()->GetNumberOfElements();
      Connections = new XdmfInt64[ Length ];
      xdmfGrid->GetTopology()->GetConnectivity()->GetValues(
        0, 
        Connections, 
        Length);
      
      NumberOfElements = 
        xdmfGrid->GetTopology()->GetShapeDesc()->GetNumberOfElements();
      ctp = cell_types = new int[ NumberOfElements ];
      
      /* Create Cell Array */
      verts = vtkCellArray::New();
      
      /* Get the pointer */
      connections = verts->WritePointer(
        NumberOfElements,
        NumberOfElements * ( 1 + NodesPerElement ));
      
      /* Connections : N p1 p2 ... pN */
      /* i.e. Triangles : 3 0 1 2    3 3 4 5   3 6 7 8 */
      index = 0;
      for( j = 0 ; j < NumberOfElements; j++ )
Berk Geveci's avatar
Berk Geveci committed
740
        {
741
742
743
744
745
746
        *ctp++ = vType;
        *connections++ = NodesPerElement;
        for( i = 0 ; i < NodesPerElement; i++ )
          {
          *connections++ = Connections[index++];
          }
Berk Geveci's avatar
Berk Geveci committed
747
        }
748
749
750
751
752
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
      } 
    else 
      {
      // Mixed Topology
      /* Create Cell Type Array */
      vtkIdTypeArray *IdArray;
      vtkIdType RealSize;
      
      Length = xdmfGrid->GetTopology()->GetConnectivity()->GetNumberOfElements();
      Connections = new XdmfInt64[ Length ];
      xdmfGrid->GetTopology()->GetConnectivity()->GetValues(
        0, 
        Connections, 
        Length );
      NumberOfElements = 
        xdmfGrid->GetTopology()->GetShapeDesc()->GetNumberOfElements();
      ctp = cell_types = new int[ NumberOfElements ];
      
      /* Create Cell Array */
      verts = vtkCellArray::New();
      
      /* Get the pointer. Make it Big enough ... too big for now */
      // cout << "::::: Length = " << Length << endl;
      connections = verts->WritePointer(
        NumberOfElements,
        Length);
      //   Length);
      /* Connections : N p1 p2 ... pN */
      /* i.e. Triangles : 3 0 1 2    3 3 4 5   3 6 7 8 */
      index = 0;
      for( j = 0 ; j < NumberOfElements; j++ )
Berk Geveci's avatar
Berk Geveci committed
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
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
        switch ( Connections[index++] )
          {
          case  XDMF_POLYVERTEX :
            vType = VTK_POLY_VERTEX;
            NodesPerElement = Connections[index++];
            break;
          case  XDMF_POLYLINE :
            vType = VTK_POLY_LINE;
            NodesPerElement = Connections[index++];
            break;
          case  XDMF_POLYGON :
            vType = VTK_POLYGON;
            NodesPerElement = Connections[index++];
            break;
          case  XDMF_TRI :
            vType = VTK_TRIANGLE;
            NodesPerElement = 3;
            break;
          case  XDMF_QUAD :
            vType = VTK_QUAD;
            NodesPerElement = 4;
            break;
          case  XDMF_TET :
            vType = VTK_TETRA;
            NodesPerElement = 4;
            break;
          case  XDMF_PYRAMID :
            vType = VTK_PYRAMID;
            NodesPerElement = 5;
            break;
          case  XDMF_WEDGE :
            vType = VTK_WEDGE;
            NodesPerElement = 6;
            break;
          case  XDMF_HEX :
            vType = VTK_HEXAHEDRON;
            NodesPerElement = 8;
            break;
          case  XDMF_EDGE_3 :
            vType = VTK_QUADRATIC_EDGE ;
            NodesPerElement = 3;
            break;
          case  XDMF_TRI_6 :
            vType = VTK_QUADRATIC_TRIANGLE ;
            NodesPerElement = 6;
            break;
          case  XDMF_QUAD_8 :
            vType = VTK_QUADRATIC_QUAD ;
            NodesPerElement = 8;
            break;
          case  XDMF_TET_10 :
            vType = VTK_QUADRATIC_TETRA ;
            NodesPerElement = 10;
            break;
          case  XDMF_PYRAMID_13 :
            vType = VTK_QUADRATIC_PYRAMID ;
            NodesPerElement = 13;
            break;
          case  XDMF_WEDGE_15 :
            vType = VTK_QUADRATIC_WEDGE ;
            NodesPerElement = 15;
            break;
          case  XDMF_HEX_20 :
            vType = VTK_QUADRATIC_HEXAHEDRON ;
            NodesPerElement = 20;
            break;
          default :
            XdmfErrorMessage("Unknown Topology Type");
            return 1;
          }
        *ctp++ = vType;
        *connections++ = NodesPerElement;
        for( i = 0 ; i < NodesPerElement; i++ )
          {
          *connections++ = Connections[index++];
          }
Berk Geveci's avatar
Berk Geveci committed
856
        }
857
858
859
860
861
862
      // Resize the Array to the Proper Size
      IdArray = verts->GetData();
      RealSize = index - 1;
      vtkDebugWithObjectMacro(this->Reader, 
                              "Resizing to " << RealSize << " elements");
      IdArray->Resize(RealSize);
Berk Geveci's avatar
Berk Geveci committed
863
      }
864

865
866
867
868
869
870
871
    delete [] Connections;
    vGrid->SetCells(cell_types, verts);
    /* OK, because of reference counting */
    verts->Delete();
    delete [] cell_types;
    vGrid->Modified();
    }  // if( xdmfGrid->GetClass() == XDMF_UNSTRUCTURED ) 
Jerry Clarke's avatar
Jerry Clarke committed
872
873
  else if( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DSMESH ||
           xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DSMESH )
874
    {
Berk Geveci's avatar
Berk Geveci committed
875
876
    vtkDebugWithObjectMacro(this->Reader, 
                            << "Setting Extents for vtkStructuredGrid");
877
878
879
    vtkStructuredGrid  *vGrid = vtkStructuredGrid::SafeDownCast(output);
    vGrid->SetExtent(upext);    
    } 
Jerry Clarke's avatar
Jerry Clarke committed
880
881
  else if (xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DCORECTMESH ||
           xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DCORECTMESH ) 
882
883
884
885
886
887
888
889
    {
#ifdef USE_IMAGE_DATA
    vtkImageData *idata = vtkImageData::SafeDownCast(output);
#else
    vtkUniformGrid *idata = vtkUniformGrid::SafeDownCast(output);
#endif
    idata->SetExtent(upext);
    }
Jerry Clarke's avatar
Jerry Clarke committed
890
891
  else if ( xdmfGrid->GetTopology()->GetTopologyType() == XDMF_2DRECTMESH ||
            xdmfGrid->GetTopology()->GetTopologyType() == XDMF_3DRECTMESH )
892
893
894
895
896
897
898
    {
    vtkRectilinearGrid *vGrid = vtkRectilinearGrid::SafeDownCast(output);
    vGrid->SetExtent(upext);    
    }
  else
    {
    vtkErrorWithObjectMacro(this->Reader,"Do not understand topology type: " 
Jerry Clarke's avatar
Jerry Clarke committed
899
                            << xdmfGrid->GetTopology()->GetTopologyType());
900
901
902
    }
  // Read Geometry
  if( ( Geometry->GetGeometryType() == XDMF_GEOMETRY_X_Y_Z ) ||
903
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_XYZ ) ||
904
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_X_Y ) ||
905
      ( Geometry->GetGeometryType() == XDMF_GEOMETRY_XY ) )
906
907
908
909
910
911
912
913
    {
    XdmfInt64   Length;
    vtkPoints   *Points;
    vtkPointSet *Pointset = vtkPointSet::SafeDownCast(output);
    
    // Special flag, for structured data
    int structured_grid = 0;
    if ( vtkStructuredGrid::SafeDownCast(output) )
914
      {
915
916
      structured_grid = 1;
      }
917
    
918
919
920
921
922
923
924
925
926
    Points = Pointset->GetPoints();
    if( !Points )
      {
      vtkDebugWithObjectMacro(this->Reader,<<  "Creating vtkPoints" );
      Points = vtkPoints::New();
      Pointset->SetPoints( Points );
      // OK Because of Reference Counting
      Points->Delete();
      }
927
    
928
929
930
    if( Geometry->GetPoints() )
      {
      if( Points )
931
        {
932
        if ( Geometry->GetPoints()->GetNumberType() == XDMF_FLOAT32_TYPE )
933
          {
934
          if ( Points->GetData()->GetDataType() != VTK_FLOAT)
935
            {
936
937
938
939
            vtkFloatArray* da = vtkFloatArray::New();
            da->SetNumberOfComponents(3);
            Points->SetData(da);
            da->Delete();
940
            }
941
942
943
944
          }
        else // means == XDMF_FLOAT64_TYPE
          {
          if ( Points->GetData()->GetDataType() != VTK_DOUBLE )
945
            {
946
947
948
949
            vtkDoubleArray* da = vtkDoubleArray::New();
            da->SetNumberOfComponents(3);
            Points->SetData(da);
            da->Delete();
950
            }
951
952
953
          }
        
        Length = Geometry->GetPoints()->GetNumberOfElements();
Berk Geveci's avatar
Berk Geveci committed
954
955
956
957
        vtkDebugWithObjectMacro(this->Reader, 
                                << "Setting Array of " << (int)Length << " = " 
                                << (int)Geometry->GetNumberOfPoints() 
                                << " Points");
958
959
960
961
962
963
964
965
966
967
        vtkIdType iskip[3] = { 0, 0, 0 };
        vtkIdType eskip[3] = { 0, 0, 0 };
        int strides_or_extents = 0;
        if ( structured_grid )
          {
          XdmfInt64     ii, jj, kk;
          vtkIdType numpoints = Geometry->GetNumberOfPoints();
          vtkIdType newnumpoints = ((upext[5] - upext[4] + 1) * (upext[3] - upext[2] + 1) * (upext[1] - upext[0] + 1));
          int cnt = 0;
          for (kk = upext[4]; kk <= upext[5]; kk ++ )
968
            {
969
            for ( jj = upext[2]; jj <= upext[3]; jj ++ )
970
              {
971
              for ( ii = upext[0]; ii <= upext[1]; ii ++ )
972
                {
973
                cnt ++;
974
975
976
                }
              }
            }
977
978
979
980
981
982
983
984
985
986
987
988
989
990
          newnumpoints = cnt;
          
          Points->SetNumberOfPoints(newnumpoints);
          vtkIdType dims[3];
          dims[0] = whext[1] - whext[0] + 1;
          dims[1] = whext[3] - whext[2] + 1;
          dims[2] = whext[5] - whext[4] + 1;
          iskip[0] = upext[0];
          iskip[1] = upext[2] * dims[0];
          iskip[2] = upext[4] * dims[0] * dims[1];
          eskip[0] = whext[1] - upext[1];
          eskip[1] = (whext[3] - upext[3]) * dims[0];
          eskip[2] = (whext[5] - upext[5]) * dims[0] * dims[1];
          if ( newnumpoints != numpoints )
991
            {
992
            strides_or_extents = 1;
993
            }
994
995
996
997
998
999
1000
          }
        else
          {
          // Unstructured grid
          Points->SetNumberOfPoints( Geometry->GetNumberOfPoints() );
          int kk;
          for ( kk = 0; kk < 6; kk ++ )
For faster browsing, not all history is shown. View entire blame