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);
Dave Demarle's avatar
Dave Demarle committed
87
vtkCxxRevisionMacro(vtkXdmfReader, "1.25");
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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  ~vtkXdmfReaderActualGrid() 
  {
  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
230
231
232
233
  int Enabled;
  vtkXdmfReaderGrid* Grid;
  vtkXdmfReaderGridCollection* Collection;
};

234
//============================================================================
Andy Cedilnik's avatar
Andy Cedilnik committed
235
class vtkXdmfReaderInternal
236
237
{
public:
238
239
  vtkXdmfReaderInternal()
    {
Berk Geveci's avatar
Berk Geveci committed
240
241
242
      this->DataItem = 0;
      this->ArrayConverter = vtkXdmfDataArray::New();
      this->NumberOfEnabledActualGrids=0;
243
244
245
    }
  ~vtkXdmfReaderInternal()
    {
Berk Geveci's avatar
Berk Geveci committed
246
247
248
249
250
251
252
253
254
      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;
255
256
    }

257
  typedef vtkstd::vector<vtkstd::string> StringListType;
258
259
  StringListType DomainList;
  XdmfXmlNode DomainPtr;
260
261

  int RequestSingleGridData(const char* currentGridName,
262
263
                            vtkXdmfReaderGrid *grid,
                            vtkInformation *outInfo,
264
265
                            vtkDataObject *output,
                            int isSubBlock);
266
267
268
269
  
  int RequestActualGridData(const char* currentGridName,
                            vtkXdmfReaderActualGrid* currentActualGrid,
                            int outputGrid,
270
                            int numberOfGrids,
271
272
273
274
275
276
277
278
                            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
279
                                   int numberOfGrids,
280
281
282
                                   vtkInformationVector* outputVector);
  
  typedef vtkstd::map<vtkstd::string,vtkXdmfReaderActualGrid> MapOfActualGrids;
283
  MapOfActualGrids ActualGrids;
284
285
  int NumberOfEnabledActualGrids;
  
286
287
288
  vtkXdmfReaderGridCollection* GetCollection(const char* collectionName);
  vtkXdmfReaderActualGrid* GetGrid(const char* gridName);
  vtkXdmfReaderActualGrid* GetGrid(int idx);
289
  vtkXdmfReaderGrid *GetXdmfGrid(const char *gridName,
290
291
                                 const char *collectionName,
                                 const char *levelName);
292
293

  vtkXdmfReader* Reader;
Jerry Clarke's avatar
Jerry Clarke committed
294
  XdmfDataItem *DataItem;
295
296
297

  // For converting arrays from XDMF to VTK format
  vtkXdmfDataArray *ArrayConverter;
298
299
};

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

309
  vtkXdmfReaderActualGrid* actualGrid = &this->ActualGrids[collectionName];
310

311
312
313
314
315
316
317
318
319
  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;
    }
320

321
322
323
324
  return actualGrid->Collection;
}

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

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

  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;
    }
373
  grid->Grid = new vtkXdmfReaderGrid; 
374
375
376
377
378
379
380
381
382
383
384
385
  return grid->Grid;
}

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

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

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

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

427
428
429
430
431
432
  unsigned int numberOfDataSets=currentActualGrid->Collection->Grids.size();
    
  currentActualGrid->Collection->UpdateCounts();
  int levels=currentActualGrid->Collection->GetNumberOfLevels();
  // mgd->SetNumberOfLevels(levels);
  // int levels = 1;
433
  
434
435
436
437
438
439
440
  int level=0;
  mgd->SetNumberOfDataSets(outputGrid, currentActualGrid->Collection->GetNumberOfDataSets(level));
  while(level<levels)
    {
    // mgd->SetNumberOfDataSets(level,currentActualGrid->Collection->GetNumberOfDataSets(level));
    ++level;
    }
441
  
442
443
444
445
446
447
448
449
450
451
  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
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
    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)
487
      {
488
489
      // mgd->SetDataSet(level,index,0); // empty, on another processor
      mgd->SetDataSet(outputGrid, index, 0); // empty, on another processor
490
      }
491
    else
492
      {
493
494
      XdmfGrid *xdmfGrid=subgrid->XMGrid;
      if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
495
        {
496
497
498
499
500
501
502
503
504
505
506
507
508
        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();
509
        }
510
511
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
      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);
538
      vtkInformation *subInfo=compInfo->GetInformation(outputGrid,index); 
539
      result=this->RequestSingleGridData("",gridIt->second,subInfo,ds,1);
540
      }
541
542
543
544
    ++currentIndex[level];
    ++gridIt;
    ++datasetIdx;
    this->Reader->UpdateProgress(1.0 * datasetIdx / numberOfDataSets);
Andy Cedilnik's avatar
Andy Cedilnik committed
545
    }
546
  return result;
Jerry Clarke's avatar
Jerry Clarke committed
547
548
}

549
//----------------------------------------------------------------------------
550
551
552
553
554
555
int vtkXdmfReaderInternal::RequestSingleGridData(
  const char* currentGridName,
  vtkXdmfReaderGrid *grid,
  vtkInformation *outInfo,
  vtkDataObject *output,
  int isSubBlock)
556
{
557
558
559
560
561
562
563
564
565
566
567
568
  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 )
569
    {
570
571
    grid->DataDescription = xdmfGrid->GetTopology()->GetShapeDesc();
    //continue;
572
573
    }
  
574
575
576
577
578
579
580
581
  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())
582
    {
583
584
    case XDMF_2DSMESH: case XDMF_2DRECTMESH: case XDMF_2DCORECTMESH:
      globalrank = 2;
585
    }
Jerry Clarke's avatar
Jerry Clarke committed
586
  if ( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED )
587
588
589
590
591
592
    {
    globalrank = 1;
    }
  
  
  vtkIdType cc;
Jerry Clarke's avatar
Jerry Clarke committed
593
594
  XdmfXmlNode attrNode;
  XdmfXmlNode dataNode;
595
596
597
598
599
600
601
602
  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
603
  if( xdmfGrid->GetTopology()->GetClass() != XDMF_UNSTRUCTURED)
604
    {
605
    if(isSubBlock)
606
      {
607
608
609
      // the composite pipeline does not set update extent, so just take
      // the whole extent for as update extent.
      outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), upext);
610
      }
611
    else
612
      {
613
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), upext);
614
      }
615
    
616
617
    outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), whext);
    
618
619
620
    start[2] = vtkMAX(0, upext[0]);
    start[1] = vtkMAX(0, upext[2]);
    start[0] = vtkMAX(0, upext[4]);
621
    
622
623
624
    count[2] = upext[1] - upext[0];
    count[1] = upext[3] - upext[2];
    count[0] = upext[5] - upext[4];
625
    }
626
  
627
  XdmfGeometry  *Geometry = xdmfGrid->GetGeometry();
628
629
  
  
630
  // Read Topology for Unstructured Grid
Jerry Clarke's avatar
Jerry Clarke committed
631
  if( xdmfGrid->GetTopology()->GetClass() == XDMF_UNSTRUCTURED ) 
632
633
634
635
636
637
638
639
640
641
642
    {
    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
643
644
645
    vtkDebugWithObjectMacro(
      this->Reader,<< "Unstructured Topology is " 
      << xdmfGrid->GetTopology()->GetTopologyTypeAsString());
Jerry Clarke's avatar
Jerry Clarke committed
646
    switch ( xdmfGrid->GetTopology()->GetTopologyType() )
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
      {
      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
668
        break;
669
670
671
672
673
674
      case  XDMF_WEDGE :
        vType = VTK_WEDGE;
        break;
      case  XDMF_HEX :
        vType = VTK_HEXAHEDRON;
        break;
Jerry Clarke's avatar
Jerry Clarke committed
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
      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
696
697
698
      case XDMF_MIXED :
        vType = -1;
        break;
699
      default :
Berk Geveci's avatar
Berk Geveci committed
700
701
        XdmfErrorMessage("Unknown Topology Type = " 
                         << xdmfGrid->GetTopology()->GetTopologyType());
702
        return 1;
703
      }
704
    if( xdmfGrid->GetTopology()->GetTopologyType() != XDMF_MIXED)
Berk Geveci's avatar
Berk Geveci committed
705
      {
706
707
      NodesPerElement = xdmfGrid->GetTopology()->GetNodesPerElement();
      if ( xdmfGrid->GetTopology()->GetConnectivity()->GetRank() == 2 )
Berk Geveci's avatar
Berk Geveci committed
708
        {
709
710
        NodesPerElement = 
          xdmfGrid->GetTopology()->GetConnectivity()->GetDimension(1);
Berk Geveci's avatar
Berk Geveci committed
711
        }
Jerry Clarke's avatar
Jerry Clarke committed
712
    
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
      /* 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
737
        {
738
739
740
741
742
743
        *ctp++ = vType;
        *connections++ = NodesPerElement;
        for( i = 0 ; i < NodesPerElement; i++ )
          {
          *connections++ = Connections[index++];
          }
Berk Geveci's avatar
Berk Geveci committed
744
        }
745
746
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
      } 
    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
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
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
        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
853
        }
854
855
856
857
858
859
      // 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
860
      }
861

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