moab
|
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