moab
DataSetConverter.h
Go to the documentation of this file.
00001 #ifndef __smoab_DataSetConverter_h
00002 #define __smoab_DataSetConverter_h
00003 
00004 #include "SimpleMoab.h"
00005 #include "CellTypeToType.h"
00006 #include "MixedCellConnectivity.h"
00007 
00008 #include <vtkCellArray.h>
00009 #include <vtkCellData.h>
00010 #include <vtkDoubleArray.h>
00011 #include <vtkFieldData.h>
00012 #include <vtkIntArray.h>
00013 #include <vtkIdTypeArray.h>
00014 #include <vtkNew.h>
00015 #include <vtkPointData.h>
00016 #include <vtkPoints.h>
00017 #include <vtkUnsignedCharArray.h>
00018 #include <vtkUnstructuredGrid.h>
00019 
00020 #include <algorithm>
00021 
00022 namespace smoab{
00023 class DataSetConverter
00024 {
00025   const smoab::Interface& Interface;
00026   moab::Interface* Moab;
00027   const smoab::Tag *Tag;
00028   bool ReadMaterialIds;
00029   bool ReadProperties;
00030   std::string MaterialName;
00031 
00032 public:
00033   DataSetConverter(const smoab::Interface& interface, const smoab::Tag* tag):
00034     Interface(interface),
00035     Moab(interface.Moab),
00036     Tag(tag),
00037     ReadMaterialIds(false),
00038     ReadProperties(false),
00039     MaterialName("Material")
00040     {
00041     }
00042 
00043   void readMaterialIds(bool add) { this->ReadMaterialIds = add; }
00044   bool readMaterialIds() const { return this->ReadMaterialIds; }
00045 
00046   void materialIdName(const std::string& name) { this->MaterialName = name; }
00047   const std::string& materialIdName() const { return this->MaterialName; }
00048 
00049   void readProperties(bool readProps) { this->ReadProperties = readProps; }
00050   bool readProperties() const { return this->ReadProperties; }
00051 
00052   //----------------------------------------------------------------------------
00053   //given a range of entity handles merge them all into a single unstructured
00054   //grid. Currently doesn't support reading properties.
00055   //Will read in material ids,  if no material id is assigned to an entity,
00056   //its cells will be given an unique id
00057   bool fill(const smoab::Range& entities,
00058             vtkUnstructuredGrid* grid) const
00059     {
00060     //create a helper datastructure which can determines all the unique point ids
00061     //and converts moab connecitvity info to vtk connectivity
00062     moab::Range cells;
00063 
00064     //append all the entities cells together into a single range
00065     typedef smoab::Range::const_iterator iterator;
00066     for(iterator i=entities.begin(); i!= entities.end(); ++i)
00067       {
00068       if(this->Tag->isComparable())
00069         {
00070         //if we are comparable only find the cells that match our tags dimension
00071         smoab::Range entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());
00072         cells.insert(entitiesCells.begin(),entitiesCells.end());
00073         }
00074       else
00075         {
00076         //this is a bad representation of all other tags, but we are presuming that
00077         //neuman and dirichlet are on entitysets with no children
00078         this->Moab->get_entities_by_handle(*i,cells);
00079         }
00080       }
00081 
00082     smoab::Range points;
00083     this->loadCellsAndPoints(cells,points,grid);
00084 
00085     if(this->readMaterialIds())
00086       {
00087       typedef std::vector<smoab::EntityHandle>::const_iterator EntityHandleIterator;
00088       typedef std::vector<int>::const_iterator IdConstIterator;
00089       typedef std::vector<int>::iterator IdIterator;
00090 
00091       std::vector<smoab::EntityHandle> searchableCells;
00092       searchableCells.reserve(grid->GetNumberOfCells());
00093       std::copy(cells.begin(),cells.end(),std::back_inserter(searchableCells));
00094       cells.clear(); //release memory we don't need
00095 
00096 
00097       std::vector<int> materialIds(entities.size());
00098       //first off iterate the entities and determine which ones
00099       //have moab material ids
00100 
00101         //wrap this area with scope, to remove local variables
00102       {
00103         smoab::MaterialTag tag;
00104         IdIterator materialIndex = materialIds.begin();
00105         for(iterator i=entities.begin();
00106             i != entities.end();
00107             ++i, ++materialIndex)
00108           {
00109           moab::Tag mtag = this->Interface.getMoabTag(tag);
00110 
00111           int value=-1;
00112           this->Moab->tag_get_data(mtag,&(*i),1,&value);
00113           *materialIndex=static_cast<int>(value);
00114           }
00115 
00116         //now determine ids for all entities that don't have materials
00117         IdConstIterator maxPos = std::max_element(materialIds.begin(),
00118                                                  materialIds.end());
00119         int maxMaterial = *maxPos;
00120         for(IdIterator i=materialIds.begin(); i!= materialIds.end(); ++i)
00121           {
00122           if(*i==-1)
00123             {
00124             *i = ++maxMaterial;
00125             }
00126           }
00127       }
00128 
00129       //now we create the material field, and set all the values
00130       vtkNew<vtkIntArray> materialSet;
00131       materialSet->SetName(this->materialIdName().c_str());
00132       materialSet->SetNumberOfValues(grid->GetNumberOfCells());
00133 
00134       IdConstIterator materialValue = materialIds.begin();
00135       for(iterator i=entities.begin(); i!= entities.end(); ++i, ++materialValue)
00136         {
00137         //this is a time vs memory trade off, I don't want to store
00138         //the all the cell ids twice over, lets use more time
00139         smoab::Range entitiesCells;
00140         if(this->Tag->isComparable())
00141           {entitiesCells = this->Interface.findEntitiesWithDimension(*i,Tag->value());}
00142         else
00143           {this->Moab->get_entities_by_handle(*i,entitiesCells);}
00144 
00145         EntityHandleIterator s_begin = searchableCells.begin();
00146         EntityHandleIterator s_end = searchableCells.end();
00147         for(iterator j=entitiesCells.begin(); j != entitiesCells.end();++j)
00148           {
00149           EntityHandleIterator result = std::lower_bound(s_begin,
00150                                                          s_end,
00151                                                          *j);
00152           std::size_t newId = std::distance(s_begin,result);
00153           materialSet->SetValue(static_cast<int>(newId), *materialValue);
00154           }
00155         }
00156 
00157       grid->GetCellData()->AddArray(materialSet.GetPointer());
00158 
00159       }
00160 
00161     return true;
00162     }
00163 
00164   //----------------------------------------------------------------------------
00165   //given a single entity handle create a unstructured grid from it.
00166   //optional third parameter is the material id to use if readMaterialIds
00167   //is on, and no material sparse tag is found for this entity
00168   bool fill(const smoab::EntityHandle& entity,
00169             vtkUnstructuredGrid* grid,
00170             const int materialId=0) const
00171     {
00172     //create a helper datastructure which can determines all the unique point ids
00173     //and converts moab connecitvity info to vtk connectivity
00174     moab::Range cells;
00175     if(this->Tag->isComparable())
00176       {
00177       //if we are comparable only find the cells that match our tags dimension
00178       cells = this->Interface.findEntitiesWithDimension(entity,Tag->value());
00179       }
00180     else
00181       {
00182       //this is a bad representation of all other tags, but we are presuming that
00183       //neuman and dirichlet are on entitysets with no children
00184       this->Moab->get_entities_by_handle(entity,cells);
00185       }
00186 
00187     smoab::Range points;
00188     this->loadCellsAndPoints(cells,points,grid);
00189 
00190 
00191     if(this->readProperties())
00192       {
00193       this->readProperties(cells,grid->GetCellData());
00194       this->readProperties(points,grid->GetPointData());
00195       }
00196 
00197     if(this->readMaterialIds())
00198       {
00199       this->readSparseTag(smoab::MaterialTag(),entity,
00200                           grid->GetNumberOfCells(),
00201                           grid->GetCellData(),
00202                           materialId);
00203       }
00204     return true;
00205     }
00206 
00207   //----------------------------------------------------------------------------
00208   void loadCellsAndPoints(const smoab::Range& cells,
00209                           smoab::Range& points,
00210                           vtkUnstructuredGrid* grid) const
00211     {
00212 
00213     smoab::MixedCellConnectivity mixConn(cells,this->Moab);
00214 
00215     //now that mixConn has all the cells properly stored, lets fixup
00216     //the ids so that they start at zero and keep the same logical ordering
00217     //as before.
00218     vtkIdType numCells, connLen;
00219     mixConn.compactIds(numCells,connLen);
00220     this->setGridsTopology(mixConn,grid,numCells,connLen);
00221 
00222     mixConn.moabPoints(points);
00223 
00224     vtkNew<vtkPoints> newPoints;
00225     this->addCoordinates(points,newPoints.GetPointer());
00226     grid->SetPoints(newPoints.GetPointer());
00227 
00228     }
00229 
00230   //----------------------------------------------------------------------------
00231   void addCoordinates(smoab::Range pointEntities, vtkPoints* pointContainer) const
00232     {
00233     //since the smoab::range are always unique and sorted
00234     //we can use the more efficient coords_iterate
00235     //call in moab, which returns moab internal allocated memory
00236     pointContainer->SetDataTypeToDouble();
00237     pointContainer->SetNumberOfPoints(pointEntities.size());
00238 
00239     //need a pointer to the allocated vtkPoints memory so that we
00240     //don't need to use an extra copy and we can bypass all vtk's check
00241     //on out of bounds
00242     double *rawPoints = static_cast<double*>(pointContainer->GetVoidPointer(0));
00243     this->Moab->get_coords(pointEntities,rawPoints);
00244     }
00245 
00246 private:
00247   //----------------------------------------------------------------------------
00248   void readProperties(smoab::Range const& entities,
00249                            vtkFieldData* field) const
00250     {
00251     if(entities.empty()) { return; }
00252 
00253     //so we get all the tags and parse out the sparse and dense tags
00254     //that we support
00255     typedef std::vector<moab::Tag>::const_iterator iterator;
00256     std::vector<moab::Tag> tags;
00257     this->Moab->tag_get_tags_on_entity(entities.front(),tags);
00258 
00259     this->readDenseTags(tags,entities,field);
00260     }
00261 
00262   //----------------------------------------------------------------------------
00263   bool readSparseTag(smoab::Tag tag,
00264                       smoab::EntityHandle const& entity,
00265                       vtkIdType length,
00266                       vtkFieldData* field,
00267                       vtkIdType defaultValue) const
00268     {
00269 
00270     typedef std::vector<moab::Tag>::const_iterator iterator;
00271     moab::Tag mtag = this->Interface.getMoabTag(tag);
00272 
00273     int value=0;
00274     moab::ErrorCode rval = this->Moab->tag_get_data(mtag,&entity,1,&value);
00275     if(rval!=moab::MB_SUCCESS)
00276       {
00277       value = defaultValue;
00278       }
00279 
00280     vtkNew<vtkIntArray> materialSet;
00281     materialSet->SetNumberOfValues(length);
00282     materialSet->SetName(this->materialIdName().c_str());
00283 
00284     int *raw = static_cast<int*>(materialSet->GetVoidPointer(0));
00285     std::fill(raw,raw+length,value);
00286 
00287     field->AddArray(materialSet.GetPointer());
00288 
00289     return true;
00290     }
00291 
00292   //----------------------------------------------------------------------------
00293   void readDenseTags(std::vector<moab::Tag> &tags,
00294                      smoab::Range const& entities,
00295                      vtkFieldData* field) const
00296     {
00297     typedef std::vector<moab::Tag>::const_iterator iterator;
00298 
00299     for(iterator i=tags.begin();i!=tags.end();++i)
00300       {
00301       moab::TagType tagType;
00302       moab::DataType tagDataType;
00303 
00304       this->Moab->tag_get_type(*i,tagType);
00305       this->Moab->tag_get_data_type(*i,tagDataType);
00306 
00307       //make sure it is only dense
00308       if(tagType != moab::MB_TAG_DENSE)
00309         {
00310         continue;
00311         }
00312       //and only integer and double
00313       if(tagDataType != moab::MB_TYPE_DOUBLE &&
00314          tagDataType != moab::MB_TYPE_INTEGER)
00315         {
00316         //unsupported type, skip to next tag
00317         continue;
00318         }
00319 
00320       //read the name of the tag
00321       std::string name;
00322       name.reserve(32);
00323       this->Moab->tag_get_name(*i,name);
00324 
00325       //read the number of components of the tag
00326       int numComps = 1;
00327 
00328       this->Moab->tag_get_length(*i,numComps);
00329 
00330       //read the data if it is one of the two types we support
00331       int size = entities.size();
00332       if(tagDataType == moab::MB_TYPE_DOUBLE)
00333         {
00334         vtkNew<vtkDoubleArray> array;
00335         array->SetName(name.c_str());
00336         array->SetNumberOfComponents(numComps);
00337         array->SetNumberOfTuples(size);
00338 
00339         //read directly into the double array
00340         this->Moab->tag_get_data(*i,entities,
00341                                  array->GetVoidPointer(0));
00342         field->AddArray(array.GetPointer());
00343         }
00344       else if(tagDataType == moab::MB_TYPE_INTEGER)
00345         {
00346         vtkNew<vtkIntArray> array;
00347         array->SetName(name.c_str());
00348         array->SetNumberOfComponents(numComps);
00349         array->SetNumberOfTuples(size);
00350 
00351         //read directly into the double array
00352         this->Moab->tag_get_data(*i,entities,
00353                                  array->GetVoidPointer(0));
00354         field->AddArray(array.GetPointer());
00355         }
00356       else
00357         {
00358         }
00359       }
00360     }
00361 
00362   //----------------------------------------------------------------------------
00363   void setGridsTopology(smoab::MixedCellConnectivity const& mixedCells,
00364                 vtkUnstructuredGrid* grid,
00365                 vtkIdType numCells,
00366                 vtkIdType numConnectivity) const
00367     {
00368     //correct the connectivity size to account for the vtk padding
00369     const vtkIdType vtkConnectivity = numCells + numConnectivity;
00370 
00371     vtkNew<vtkIdTypeArray> cellArray;
00372     vtkNew<vtkIdTypeArray> cellLocations;
00373     vtkNew<vtkUnsignedCharArray> cellTypes;
00374 
00375     cellArray->SetNumberOfValues(vtkConnectivity);
00376     cellLocations->SetNumberOfValues(numCells);
00377     cellTypes->SetNumberOfValues(numCells);
00378 
00379     vtkIdType* rawArray = static_cast<vtkIdType*>(cellArray->GetVoidPointer(0));
00380     vtkIdType* rawLocations = static_cast<vtkIdType*>(cellLocations->GetVoidPointer(0));
00381     unsigned char* rawTypes = static_cast<unsigned char*>(cellTypes->GetVoidPointer(0));
00382 
00383     mixedCells.copyToVtkCellInfo(rawArray,rawLocations,rawTypes);
00384 
00385     vtkNew<vtkCellArray> cells;
00386     cells->SetCells(numCells,cellArray.GetPointer());
00387     grid->SetCells(cellTypes.GetPointer(),
00388                    cellLocations.GetPointer(),
00389                    cells.GetPointer(),
00390                    NULL,NULL);
00391     }
00392 };
00393 
00394 }
00395 #endif // __smoab_DataSetConverter_h
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines