1
  2
  3
  4
  5
  6
  7
  8
  9
 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
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
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
#include "vtkMoabReader.h"

#include "SimpleMoab.h"
#include "DataSetConverter.h"


#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkUnstructuredGrid.h"

#include "vtkPolyData.h"


vtkStandardNewMacro(vtkMoabReader)
//------------------------------------------------------------------------------
vtkMoabReader::vtkMoabReader()
  {
  this->SetNumberOfInputPorts(0);
  this->FileName = NULL;
  }

//------------------------------------------------------------------------------
vtkMoabReader::~vtkMoabReader()
  {
  }

//------------------------------------------------------------------------------
int vtkMoabReader::RequestInformation(vtkInformation *request,
                       vtkInformationVector **inputVector,
                       vtkInformationVector *outputVector)
{

  //todo. Walk the file and display all the 2d and 3d elements that the users
  //could possibly want to load
  return this->Superclass::RequestInformation(request,inputVector,outputVector);
}

//------------------------------------------------------------------------------
int vtkMoabReader::RequestData(vtkInformation *vtkNotUsed(request),<--- The function 'RequestData' is never used.
                vtkInformationVector **vtkNotUsed(inputVector),
                vtkInformationVector *outputVector)
{
  //First pass is lets load in all 3d elements in a block called Volumes,
  //and load all 2d elements in a block called Surfaces

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

  vtkNew<vtkMultiBlockDataSet> volumeRoot;
  vtkNew<vtkMultiBlockDataSet> boundaryRoot;
  vtkNew<vtkMultiBlockDataSet> surfaceRoot;
  vtkNew<vtkMultiBlockDataSet> materialRoot;
  vtkNew<vtkMultiBlockDataSet> neumannRoot;
  vtkNew<vtkMultiBlockDataSet> dirichletRoot;

  const int blockIndex = output->GetNumberOfBlocks();
  output->SetBlock(blockIndex,volumeRoot.GetPointer());
  output->SetBlock(blockIndex+1,boundaryRoot.GetPointer());
  output->SetBlock(blockIndex+2,surfaceRoot.GetPointer());
  output->SetBlock(blockIndex+3,materialRoot.GetPointer());
  output->SetBlock(blockIndex+4,neumannRoot.GetPointer());
  output->SetBlock(blockIndex+5,dirichletRoot.GetPointer());

  //boring work, set the names of the blocks
  output->GetMetaData(blockIndex)->Set(vtkCompositeDataSet::NAME(), "Volumes");
  output->GetMetaData(blockIndex+1)->Set(vtkCompositeDataSet::NAME(), "Boundary");
  output->GetMetaData(blockIndex+2)->Set(vtkCompositeDataSet::NAME(), "Surfaces");
  output->GetMetaData(blockIndex+3)->Set(vtkCompositeDataSet::NAME(), "Material");
  output->GetMetaData(blockIndex+4)->Set(vtkCompositeDataSet::NAME(), "Neumann Sets");
  output->GetMetaData(blockIndex+5)->Set(vtkCompositeDataSet::NAME(), "Dirichlet Sets");

  smoab::GeomTag geom3Tag(3);
  smoab::GeomTag geom2Tag(2);
  smoab::GeomTag geom1Tag(1);
  smoab::MaterialTag matTag;
  smoab::NeumannTag neTag;
  smoab::DirichletTag diTag;

  smoab::Interface interface(this->FileName);
  this->CreateSubBlocks(volumeRoot, &interface, &geom3Tag);
  this->CreateSubBlocks(boundaryRoot, &interface, &geom3Tag, &geom2Tag);

  this->CreateSubBlocks(surfaceRoot, &interface, &geom2Tag);
  this->CreateSubBlocks(materialRoot, &interface, &matTag, &geom3Tag);
  this->CreateSubBlocks(neumannRoot, &interface, &neTag);
  this->CreateSubBlocks(dirichletRoot, &interface, &diTag);

  return 1;
}


//------------------------------------------------------------------------------
void vtkMoabReader::CreateSubBlocks(vtkNew<vtkMultiBlockDataSet> & root,
                                    smoab::Interface* interface,
                                    smoab::Tag const* parentTag,
                                    smoab::Tag const* extractTag)
{
  if(!extractTag)
    {
    extractTag = parentTag;
    }
  //basic premise: query the database for all tagged elements and create a new
  //multiblock elemenent for each
  smoab::DataSetConverter converter(*interface,extractTag);
  converter.readMaterialIds(true);
  converter.readProperties(true);

  smoab::EntityHandle rootHandle = interface->getRoot();
  smoab::Range parents = interface->findEntityRootParents(rootHandle);
  smoab::Range dimEnts = interface->findEntitiesWithTag(*parentTag,
                                                       rootHandle);

  smoab::Range geomParents = smoab::intersect(parents,dimEnts);

  parents.clear(); //remove this range as it is unneeded
  dimEnts.clear();

  //now each item in range can be extracted into a different grid
  typedef smoab::Range::iterator iterator;
  vtkIdType index = 0;
  for(iterator i=geomParents.begin(); i != geomParents.end(); ++i)
    {
    vtkNew<vtkUnstructuredGrid> block;
    //fill the dataset with geometry and properties
    converter.fill(*i, block.GetPointer(),index);

    //only add it if we have cells found
    if(block->GetNumberOfCells() > 0)
      {
      root->SetBlock(index,block.GetPointer());
      std::string name = interface->name(*i);
      if(name.size() > 0)
        {
        root->GetMetaData(index)->Set(vtkCompositeDataSet::NAME(), name.c_str());
        }
      ++index;
      }
    }
}

//------------------------------------------------------------------------------
void vtkMoabReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
}