moab
|
00001 00002 #ifndef __smoab_SimpleMoab_h 00003 #define __smoab_SimpleMoab_h 00004 00005 #include "moab/Core.hpp" 00006 #include "moab/Interface.hpp" 00007 #include "moab/Range.hpp" 00008 #include "MBTagConventions.hpp" 00009 00010 #include <iostream> 00011 00012 namespace smoab 00013 { 00014 //adjacency intersect / union named enum to match 00015 //the types from moab 00016 enum adjacency_type{INTERSECT=moab::Interface::INTERSECT, 00017 UNION=moab::Interface::UNION}; 00018 00019 //make our range equal to moabs range 00020 typedef moab::Range Range; 00021 typedef moab::EntityHandle EntityHandle; 00022 typedef moab::EntityType EntityType; 00023 00024 //bring in range functions 00025 using moab::intersect; 00026 using moab::subtract; 00027 using moab::unite; 00028 00029 //forward declare this->Moab for Tag 00030 struct Interface; 00031 00032 //forward declar the DataSetConverter so it can be a friend of Interface 00033 class DataSetConverter; 00034 00035 class Tag 00036 { 00037 const std::string Name_; 00038 public: 00039 Tag(std::string const& n):Name_(n) 00040 { 00041 } 00042 00043 virtual ~Tag() 00044 { 00045 } 00046 00047 const char* name() const { return this->Name_.c_str(); } 00048 moab::DataType virtual dataType() const { return moab::MB_TYPE_INTEGER; } 00049 virtual bool isComparable() const { return false; } 00050 virtual int value() const { return int(); } 00051 }; 00052 00053 //lightweight structs to wrap set names, so we detected 00054 //incorrect names at compile time. In the future I expect material and 00055 //boundary conditions to be comparable 00056 class MaterialTag : public Tag{ public: MaterialTag():Tag("MATERIAL_SET"){}}; 00057 class DirichletTag : public Tag{ public: DirichletTag():Tag("DIRICHLET_SET"){}}; 00058 class NeumannTag: public Tag{public: NeumannTag():Tag("NEUMANN_SET"){}}; 00059 class GroupTag: public Tag{ public: GroupTag():Tag("GROUP"){}}; 00060 00061 //geom is the only comparable tag, since it can have a dimension. 00062 class GeomTag: public Tag 00063 { 00064 int dim; 00065 public: 00066 GeomTag(int d):Tag("GEOM_DIMENSION"),dim(d){} 00067 GeomTag():Tag("GEOM_DIMENSION"), dim(0){} 00068 00069 bool isComparable() const { return dim > 0; } 00070 int value() const { return dim; } 00071 }; 00072 00073 //light weight wrapper on a moab this->Moab that exposes only the reduced class 00074 //that we need 00075 class Interface 00076 { 00077 public: 00078 Interface(const std::string &file) 00079 { 00080 this->Moab = new moab::Core(); 00081 this->Moab->load_file(file.c_str()); 00082 } 00083 00084 ~Interface() 00085 { 00086 if(this->Moab) 00087 { 00088 delete this->Moab; 00089 this->Moab = NULL; 00090 } 00091 } 00092 00093 //---------------------------------------------------------------------------- 00094 moab::Tag getMoabTag(const smoab::Tag& simpleTag) const 00095 { 00096 moab::Tag tag; 00097 this->Moab->tag_get_handle(simpleTag.name(), 00098 1, 00099 simpleTag.dataType(), 00100 tag); 00101 return tag; 00102 } 00103 00104 //---------------------------------------------------------------------------- 00105 //returns the moab name for the given entity handle if it has a sparse Name tag 00106 std::string name(const smoab::EntityHandle& entity) const 00107 { 00108 moab::Tag nameTag; 00109 moab::ErrorCode rval = this->Moab->tag_get_handle(NAME_TAG_NAME, 00110 NAME_TAG_SIZE, 00111 moab::MB_TYPE_OPAQUE, 00112 nameTag); 00113 if(rval != moab::MB_SUCCESS) { return std::string(); } 00114 00115 char name[NAME_TAG_SIZE]; 00116 rval = this->Moab->tag_get_data(nameTag,&entity,1,&name); 00117 if(rval != moab::MB_SUCCESS) { return std::string(); } 00118 00119 return std::string(name); 00120 } 00121 00122 00123 //---------------------------------------------------------------------------- 00124 smoab::EntityHandle getRoot() const { return this->Moab->get_root_set(); } 00125 00126 //---------------------------------------------------------------------------- 00127 smoab::Range findEntities(const smoab::EntityHandle root, moab::EntityType type) const 00128 { 00129 smoab::Range result; 00130 // get all the sets of that type in the mesh 00131 this->Moab->get_entities_by_type(root, type, result); 00132 return result; 00133 } 00134 00135 //---------------------------------------------------------------------------- 00136 //Find all entities with a given tag. We don't use geom as a tag as that 00137 //isn't a fast operation. Yes finding the intersection of geom entities and 00138 //a material / boundary tag will be more work, but it is rarely done currently 00139 //Returns the found group of entities 00140 smoab::Range findEntitiesWithTag (const smoab::Tag& tag, smoab::EntityHandle root, 00141 moab::EntityType type = moab::MBENTITYSET) const 00142 { 00143 smoab::Range result; 00144 00145 moab::Tag t = this->getMoabTag(tag); 00146 00147 // get all the entities of that type in the mesh 00148 this->Moab->get_entities_by_type_and_tag(root, type, &t, NULL, 1,result); 00149 00150 00151 if(tag.isComparable()) 00152 { 00153 int value=0; 00154 //now we have to remove any that doesn't match the tag value 00155 smoab::Range resultMatchingTag; 00156 typedef moab::Range::const_iterator iterator; 00157 for(iterator i=result.begin(); 00158 i != result.end(); 00159 ++i) 00160 { 00161 value = 0; 00162 moab::EntityHandle handle = *i; 00163 this->Moab->tag_get_data(t, &handle, 1, &value); 00164 if(value == tag.value()) 00165 { 00166 resultMatchingTag.insert(*i); 00167 } 00168 } 00169 00170 return resultMatchingTag; 00171 } 00172 else 00173 { 00174 //we return all the items we found 00175 return result; 00176 } 00177 } 00178 00179 //---------------------------------------------------------------------------- 00180 //Find all entities from a given root of a given dimensionality 00181 smoab::Range findEntitiesWithDimension(const smoab::EntityHandle root, 00182 int dimension) const 00183 { 00184 typedef smoab::Range::const_iterator iterator; 00185 00186 smoab::Range result; 00187 this->Moab->get_entities_by_dimension(root,dimension,result); 00188 00189 00190 smoab::Range children; 00191 this->Moab->get_child_meshsets(root,children,0); 00192 for(iterator i=children.begin(); i !=children.end();++i) 00193 { 00194 this->Moab->get_entities_by_dimension(*i,dimension,result); 00195 } 00196 return result; 00197 } 00198 00199 //---------------------------------------------------------------------------- 00200 smoab::Range findAdjacentEntities(const smoab::EntityHandle& entity, 00201 int dimension) const 00202 { 00203 const int adjType = static_cast<int>(smoab::INTERSECT); 00204 smoab::Range result; 00205 const bool create_if_missing = false; 00206 this->Moab->get_adjacencies(&entity, 00207 1, 00208 dimension, 00209 create_if_missing, 00210 result, 00211 adjType); 00212 return result; 00213 } 00214 00215 //---------------------------------------------------------------------------- 00216 smoab::Range findAdjacentEntities(const smoab::Range& range, 00217 int dimension, 00218 const smoab::adjacency_type type = smoab::UNION) const 00219 { 00220 //the smoab and moab adjacent intersection enums are in the same order 00221 const int adjType = static_cast<int>(type); 00222 smoab::Range result; 00223 const bool create_if_missing = false; 00224 this->Moab->get_adjacencies(range,dimension, 00225 create_if_missing, 00226 result, 00227 adjType); 00228 00229 return result; 00230 } 00231 00232 //---------------------------------------------------------------------------- 00233 //Find all elements in the database that have children and zero parents. 00234 //this doesn't find 00235 smoab::Range findEntityRootParents(smoab::EntityHandle const& root) const 00236 { 00237 smoab::Range parents; 00238 00239 typedef moab::Range::const_iterator iterator; 00240 moab::Range sets; 00241 00242 this->Moab->get_entities_by_type(root, moab::MBENTITYSET, sets); 00243 for(iterator i=sets.begin(); i!=sets.end();++i) 00244 { 00245 int numParents=0,numChildren=0; 00246 this->Moab->num_parent_meshsets(*i,&numParents); 00247 if(numParents==0) 00248 { 00249 this->Moab->num_child_meshsets(*i,&numChildren); 00250 if(numChildren>=0) 00251 { 00252 parents.insert(*i); 00253 } 00254 } 00255 } 00256 return parents; 00257 } 00258 00259 //---------------------------------------------------------------------------- 00260 //finds entities that have zero children and zero parents 00261 smoab::Range findDetachedEntities(moab::EntityHandle const& root) const 00262 { 00263 smoab::Range detached; 00264 00265 typedef moab::Range::const_iterator iterator; 00266 moab::Range sets; 00267 00268 this->Moab->get_entities_by_type(root, moab::MBENTITYSET, sets); 00269 for(iterator i=sets.begin(); i!=sets.end();++i) 00270 { 00271 int numParents=0,numChildren=0; 00272 this->Moab->num_parent_meshsets(*i,&numParents); 00273 if(numParents==0) 00274 { 00275 this->Moab->num_child_meshsets(*i,&numChildren); 00276 if(numChildren==0) 00277 { 00278 detached.insert(*i); 00279 } 00280 } 00281 } 00282 return detached; 00283 } 00284 00285 //---------------------------------------------------------------------------- 00286 //find all children of the entity passed in that has multiple parents 00287 smoab::Range findEntitiesWithMultipleParents(smoab::EntityHandle const& root) 00288 { 00289 smoab::Range multipleParents; 00290 typedef moab::Range::const_iterator iterator; 00291 00292 //for all the elements in the range, find all items with multiple parents 00293 moab::Range children; 00294 this->Moab->get_child_meshsets(root,children,0); 00295 for(iterator i=children.begin(); i!=children.end();++i) 00296 { 00297 int numParents=0; 00298 this->Moab->num_parent_meshsets(*i,&numParents); 00299 if(numParents>1) 00300 { 00301 multipleParents.insert(*i); 00302 } 00303 } 00304 return multipleParents; 00305 } 00306 00307 //---------------------------------------------------------------------------- 00308 //prints all elements in a range objects 00309 void printRange(smoab::Range const& range) 00310 { 00311 typedef Range::const_iterator iterator; 00312 for(iterator i=range.begin(); i!=range.end(); ++i) 00313 { 00314 std::cout << "entity id: " << *i << std::endl; 00315 this->Moab->list_entity(*i); 00316 } 00317 } 00318 00319 friend class smoab::DataSetConverter; 00320 private: 00321 moab::Interface* Moab; 00322 }; 00323 00324 } 00325 00326 #endif