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