MOAB: Mesh Oriented datABase  (version 5.2.1)
SimpleMoab.h
Go to the documentation of this file.
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 "moab/CN.hpp"
00009 
00010 #include "MBTagConventions.hpp"
00011 
00012 #include <iostream>
00013 
00014 namespace smoab
00015 {
00016 //adjacency intersect / union named enum to match
00017 //the types from moab
00018 enum adjacency_type{INTERSECT=moab::Interface::INTERSECT,
00019                     UNION=moab::Interface::UNION};
00020 
00021 //make our range equal to moabs range
00022 typedef moab::Range Range;
00023 typedef moab::EntityHandle EntityHandle;
00024 typedef moab::EntityType EntityType;
00025 
00026 //bring in range functions
00027 using moab::intersect;
00028 using moab::subtract;
00029 using moab::unite;
00030 
00031 class Tag
00032 {
00033   const std::string Name_;
00034 public:
00035   Tag(std::string const& n):Name_(n)
00036     {
00037     }
00038 
00039   virtual ~Tag()
00040     {
00041     }
00042 
00043   const char* name() const { return this->Name_.c_str(); }
00044   moab::DataType virtual dataType() const { return moab::MB_TYPE_INTEGER; }
00045   virtual bool isComparable() const { return false; }
00046   virtual int value() const { return int(); }
00047 };
00048 
00049 //lightweight structs to wrap set names, so we detected
00050 //incorrect names at compile time. In the future I expect material and
00051 //boundary conditions to be comparable
00052 class MaterialTag : public Tag{ public: MaterialTag():Tag("MATERIAL_SET"){}};
00053 class DirichletTag : public Tag{ public: DirichletTag():Tag("DIRICHLET_SET"){}};
00054 class NeumannTag: public Tag{public:  NeumannTag():Tag("NEUMANN_SET"){}};
00055 class GroupTag: public Tag{ public: GroupTag():Tag("GROUP"){}};
00056 
00057 //geom is the only comparable tag, since it can have a dimension.
00058 class GeomTag: public Tag
00059   {
00060   int dim;
00061 public:
00062   GeomTag(int d):Tag("GEOM_DIMENSION"),dim(d){}
00063   GeomTag():Tag("GEOM_DIMENSION"), dim(0){}
00064 
00065   virtual ~GeomTag(){}
00066 
00067   bool isComparable() const { return dim > 0; }
00068   int value() const { return dim; }
00069   };
00070 
00071 
00072 //forward declare this->Moab for Tag
00073 struct Interface;
00074 
00075 //forward declare the DataSetConverter so it can be a friend of Interface
00076 class DataSetConverter;
00077 
00078 //forward declare the LoadGeometry so it can be a friend of Interface
00079 namespace detail{ class LoadGeometry; }
00080 namespace detail{ class LoadPoly; }
00081 
00082 
00083 //light weight wrapper on a moab this->Moab that exposes only the reduced class
00084 //that we need
00085 class Interface
00086 {
00087 public:
00088   Interface(const std::string &file)
00089     {
00090     this->Moab = new moab::Core();
00091     this->Moab->load_file(file.c_str());
00092     }
00093 
00094   ~Interface()
00095     {
00096     if(this->Moab)
00097       {
00098       delete this->Moab;
00099       this->Moab = NULL;
00100       }
00101     }
00102 
00103   //----------------------------------------------------------------------------
00104   moab::Tag getMoabTag(const smoab::Tag& simpleTag) const
00105     {
00106     moab::Tag tag;
00107     this->Moab->tag_get_handle(simpleTag.name(),
00108                                1,
00109                                simpleTag.dataType(),
00110                                tag);
00111     return tag;
00112     }
00113 
00114   //----------------------------------------------------------------------------
00115   template<typename T>
00116   T getDefaultTagVaue(moab::Tag tag) const
00117     {
00118     T defaultValue;
00119     this->Moab->tag_get_default_value(tag,&defaultValue);
00120     return defaultValue;
00121     }
00122 
00123   //----------------------------------------------------------------------------
00124   template<typename T>
00125   T getDefaultTagVaue(smoab::Tag tag) const
00126     {
00127     return this->getDefaultTagVaue<T>(getMoabTag(tag));
00128     }
00129 
00130   //----------------------------------------------------------------------------
00131   template<typename T>
00132   T getTagData(moab::Tag tag, const smoab::EntityHandle& entity, T value) const
00133     {
00134     this->Moab->tag_get_data(tag,&entity,1,&value);
00135     return value;
00136     }
00137 
00138   //----------------------------------------------------------------------------
00139   template<typename T>
00140   T getTagData(smoab::Tag tag, const smoab::EntityHandle& entity, T value = T()) const
00141     {
00142     return this->getTagData(getMoabTag(tag),entity,value);
00143     }
00144 
00145   //----------------------------------------------------------------------------
00146   //returns the moab name for the given entity handle if it has a sparse Name tag
00147   std::string name(const smoab::EntityHandle& entity) const
00148     {
00149     moab::Tag nameTag;
00150     moab::ErrorCode rval = this->Moab->tag_get_handle(NAME_TAG_NAME,
00151                                                       NAME_TAG_SIZE,
00152                                                       moab::MB_TYPE_OPAQUE,
00153                                                       nameTag);
00154     if(rval != moab::MB_SUCCESS) { return std::string(); }
00155 
00156     char name[NAME_TAG_SIZE];
00157     rval = this->Moab->tag_get_data(nameTag,&entity,1,&name);
00158     if(rval != moab::MB_SUCCESS) { return std::string(); }
00159 
00160     return std::string(name);
00161     }
00162 
00163   //----------------------------------------------------------------------------
00164   //returns the geometeric dimension of an entity.
00165   int dimension(const smoab::EntityHandle& entity) const
00166     {
00167     return this->Moab->dimension_from_handle(entity);
00168     }
00169 
00170   //----------------------------------------------------------------------------
00171   //returns the geometeric dimension of an entity.
00172   smoab::EntityType entityType(const smoab::EntityHandle& entity) const
00173     {
00174     return this->Moab->type_from_handle(entity);
00175     }
00176 
00177   //----------------------------------------------------------------------------
00178   smoab::EntityHandle getRoot() const { return this->Moab->get_root_set(); }
00179 
00180   //----------------------------------------------------------------------------
00181   smoab::Range findEntities(const smoab::EntityHandle root, moab::EntityType type) const
00182     {
00183     smoab::Range result;
00184     // get all the sets of that type in the mesh
00185     this->Moab->get_entities_by_type(root, type, result);
00186     return result;
00187     }
00188 
00189   //----------------------------------------------------------------------------
00190   //given a single entity handle find all items in that mesh set that aren't
00191   //them selves entitysets. If recurse is true we also recurse sub entitysets
00192   smoab::Range findAllMeshEntities(smoab::EntityHandle const& entity,
00193                                    bool recurse=false) const
00194     {
00195     smoab::Range result;
00196     this->Moab->get_entities_by_handle(entity,result,recurse);
00197     return result;
00198     }
00199 
00200   //----------------------------------------------------------------------------
00201   //Find all entities with a given tag. We don't use geom as a tag as that
00202   //isn't a fast operation. Yes finding the intersection of geom entities and
00203   //a material / boundary tag will be more work, but it is rarely done currently
00204   //Returns the found group of entities
00205   smoab::Range findEntitiesWithTag (const smoab::Tag& tag, smoab::EntityHandle root,
00206                                     moab::EntityType type = moab::MBENTITYSET) const
00207     {
00208     smoab::Range result;
00209 
00210     moab::Tag t = this->getMoabTag(tag);
00211 
00212     // get all the entities of that type in the mesh
00213     this->Moab->get_entities_by_type_and_tag(root, type, &t, NULL, 1,result);
00214 
00215 
00216     if(tag.isComparable())
00217       {
00218       int value=0;
00219       //now we have to remove any that doesn't match the tag value
00220       smoab::Range resultMatchingTag;
00221       typedef moab::Range::const_iterator iterator;
00222       for(iterator i=result.begin();
00223           i != result.end();
00224           ++i)
00225         {
00226         value = 0;
00227         moab::EntityHandle handle = *i;
00228         this->Moab->tag_get_data(t, &handle, 1, &value);
00229         if(value == tag.value())
00230           {
00231           resultMatchingTag.insert(*i);
00232           }
00233         }
00234 
00235       return resultMatchingTag;
00236       }
00237     else
00238       {
00239       //we return all the items we found
00240       return result;
00241       }
00242     }
00243 
00244   //----------------------------------------------------------------------------
00245   //Find all entities from a given root of a given dimensionality
00246   smoab::Range findEntitiesWithDimension(const smoab::EntityHandle root,
00247                                          const int dimension,
00248                                          bool recurse=false) const
00249     {
00250     typedef smoab::Range::const_iterator iterator;
00251 
00252     smoab::Range result;
00253     this->Moab->get_entities_by_dimension(root,dimension,result,recurse);
00254 
00255     if(recurse)
00256       {
00257       smoab::Range children;
00258       this->Moab->get_child_meshsets(root,children,0);
00259       for(iterator i=children.begin(); i !=children.end();++i)
00260         {
00261         this->Moab->get_entities_by_dimension(*i,dimension,result);
00262         }
00263       }
00264     return result;
00265     }
00266 
00267 
00268   //----------------------------------------------------------------------------
00269   smoab::Range findHighestDimensionEntities(const smoab::EntityHandle& entity,
00270                                             bool recurse=false) const
00271     {
00272     //the goal is to load all entities that are not entity sets of this
00273     //node, while also subsetting by the highest dimension
00274 
00275     //lets find the entities of only the highest dimension
00276     int num_ents=0;
00277     int dim=3;
00278     while(num_ents<=0&&dim>0)
00279       {
00280       this->Moab->get_number_entities_by_dimension(entity,dim,num_ents,recurse);
00281       --dim;
00282       }
00283     ++dim; //reincrement to correct last decrement
00284     if(num_ents > 0)
00285       {
00286       //we have found entities of a given dimension
00287       return this->findEntitiesWithDimension(entity,dim,recurse);
00288       }
00289     return smoab::Range();
00290     }
00291 
00292   //----------------------------------------------------------------------------
00293   //Find all elements in the database that have children and zero parents.
00294   //this doesn't find
00295   smoab::Range findEntityRootParents(const smoab::EntityHandle& root) const
00296     {
00297     smoab::Range parents;
00298 
00299     typedef moab::Range::const_iterator iterator;
00300     moab::Range sets;
00301 
00302     this->Moab->get_entities_by_type(root, moab::MBENTITYSET, sets);
00303     for(iterator i=sets.begin(); i!=sets.end();++i)
00304       {
00305       int numParents=0,numChildren=0;
00306       this->Moab->num_parent_meshsets(*i,&numParents);
00307       if(numParents==0)
00308         {
00309         this->Moab->num_child_meshsets(*i,&numChildren);
00310         if(numChildren>=0)
00311           {
00312           parents.insert(*i);
00313           }
00314         }
00315       }
00316     return parents;
00317     }
00318 
00319   //----------------------------------------------------------------------------
00320   //finds entities that have zero children and zero parents
00321   smoab::Range findDetachedEntities(const moab::EntityHandle& root) const
00322     {
00323     smoab::Range detached;
00324 
00325     typedef moab::Range::const_iterator iterator;
00326     moab::Range sets;
00327 
00328     this->Moab->get_entities_by_type(root, moab::MBENTITYSET, sets);
00329     for(iterator i=sets.begin(); i!=sets.end();++i)
00330       {
00331       int numParents=0,numChildren=0;
00332       this->Moab->num_parent_meshsets(*i,&numParents);
00333       if(numParents==0)
00334         {
00335         this->Moab->num_child_meshsets(*i,&numChildren);
00336         if(numChildren==0)
00337           {
00338           detached.insert(*i);
00339           }
00340         }
00341       }
00342     return detached;
00343     }
00344 
00345   //----------------------------------------------------------------------------
00346   //find all children of the entity passed in that has multiple parents
00347   smoab::Range findEntitiesWithMultipleParents(const smoab::EntityHandle& root) const
00348     {
00349     smoab::Range multipleParents;
00350     typedef moab::Range::const_iterator iterator;
00351 
00352     //for all the elements in the range, find all items with multiple parents
00353     moab::Range children;
00354     this->Moab->get_child_meshsets(root,children,0);
00355     for(iterator i=children.begin(); i!=children.end();++i)
00356       {
00357       int numParents=0;
00358       this->Moab->num_parent_meshsets(*i,&numParents);
00359       if(numParents>1)
00360         {
00361         multipleParents.insert(*i);
00362         }
00363       }
00364     return multipleParents;
00365     }
00366 
00367   //----------------------------------------------------------------------------
00368   //find all entities that are adjacent to a single entity
00369   smoab::Range findAdjacencies(const smoab::EntityHandle& entity,
00370                                int dimension) const
00371     {
00372     const int adjType = static_cast<int>(smoab::INTERSECT);
00373     smoab::Range result;
00374     const bool create_if_missing = false;
00375     this->Moab->get_adjacencies(&entity,
00376                                 1,
00377                                 dimension,
00378                                 create_if_missing,
00379                                 result,
00380                                 adjType);
00381     return result;
00382     }
00383 
00384   //----------------------------------------------------------------------------
00385   smoab::Range findAdjacencies(const smoab::Range& range,
00386                                     int dimension,
00387                                     const smoab::adjacency_type type = smoab::UNION) const
00388     {
00389     //the smoab and moab adjacent intersection enums are in the same order
00390     const int adjType = static_cast<int>(type);
00391     smoab::Range result;
00392     const bool create_if_missing = false;
00393     this->Moab->get_adjacencies(range,dimension,
00394                                 create_if_missing,
00395                                 result,
00396                                 adjType);
00397 
00398     return result;
00399     }
00400 
00401   //----------------------------------------------------------------------------
00402   //create adjacencies, only works when the dimension requested is lower than
00403   //dimension of the range of entities
00404   smoab::Range createAdjacencies(const smoab::Range& range,
00405                                  int dimension,
00406                                  const smoab::adjacency_type type = smoab::UNION) const
00407     {
00408     //the smoab and moab adjacent intersection enums are in the same order
00409     const int adjType = static_cast<int>(type);
00410     smoab::Range result;
00411     const bool create_if_missing = true;
00412     this->Moab->get_adjacencies(range,dimension,
00413                                 create_if_missing,
00414                                 result,
00415                                 adjType);
00416 
00417     return result;
00418     }
00419 
00420   //----------------------------------------------------------------------------
00421   int numChildMeshSets(const smoab::EntityHandle& root) const
00422     {
00423     int numChildren;
00424     this->Moab->num_child_meshsets(root,&numChildren);
00425     return numChildren;
00426     }
00427 
00428   //----------------------------------------------------------------------------
00429   smoab::Range getChildSets(const smoab::EntityHandle& root) const
00430     {
00431     smoab::Range children;
00432     this->Moab->get_child_meshsets(root,children,0);
00433     return children;
00434     }
00435 
00436   //----------------------------------------------------------------------------
00437   //remove a collection of entities from the database
00438   void remove(smoab::Range const& toDelete) const
00439     {
00440     this->Moab->delete_entities(toDelete);
00441     }
00442 
00443   //----------------------------------------------------------------------------
00444   //a entityHandle with value zero means no side element was found
00445   smoab::EntityHandle sideElement(smoab::EntityHandle const& cell,
00446                                   int dim, int side) const
00447     {
00448     smoab::EntityHandle result(0);
00449     this->Moab->side_element(cell,dim,side,result);
00450     return result;
00451     }
00452 
00453   //----------------------------------------------------------------------------
00454   //returns all the existing side elements of a cell, elements that
00455   //are zero mean that side element doesn't exist
00456   std::vector<smoab::EntityHandle> sideElements(
00457                                     smoab::EntityHandle const& cell,
00458                                     int dim) const
00459     {
00460     const EntityType volumeCellType = this->Moab->type_from_handle(cell);
00461     const int numSides = static_cast<int>(moab::CN::NumSubEntities(
00462                                           volumeCellType, dim));
00463 
00464     std::vector<smoab::EntityHandle> result(numSides);
00465     for (int side = 0; side < numSides; ++side)
00466       {
00467       smoab::EntityHandle *sideElem = &result[side]; //get memory of vector
00468       this->Moab->side_element(cell,dim,side,*sideElem);
00469       }
00470     return result;
00471     }
00472 
00473   //----------------------------------------------------------------------------
00474   //prints all elements in a range objects
00475   void printRange(const smoab::Range& range) const
00476     {
00477     typedef Range::const_iterator iterator;
00478     for(iterator i=range.begin(); i!=range.end(); ++i)
00479       {
00480       std::cout << "entity id: " << *i << std::endl;
00481       this->Moab->list_entity(*i);
00482       }
00483     }
00484   friend class smoab::DataSetConverter;
00485   friend class smoab::detail::LoadGeometry;
00486   friend class smoab::detail::LoadPoly;
00487 private:
00488   moab::Interface* Moab;
00489 };
00490 
00491 //----------------------------------------------------------------------------
00492 void RangeToVector(const smoab::Range &range,
00493                    std::vector<smoab::EntityHandle>& vector )
00494 {
00495   vector.reserve(range.size());
00496   std::copy(range.begin(),
00497             range.end(),
00498             std::back_inserter(vector));
00499 }
00500 
00501 
00502 
00503 }
00504 
00505 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines