MOAB: Mesh Oriented datABase  (version 5.4.1)
FaceSets.h
Go to the documentation of this file.
00001 #ifndef __smoab_FaceSets_h
00002 #define __smoab_FaceSets_h
00003 
00004 #include "CellSets.h"
00005 #include <set>
00006 
00007 namespace smoab
00008 {
00009 //----------------------------------------------------------------------------
00010 class FaceCellSet : public CellSet
00011 {
00012   public:
00013     FaceCellSet( int id, smoab::EntityHandle p, const smoab::Range& cells ) : CellSet( p, cells ), ID( id ) {}
00014 
00015     int faceId() const
00016     {
00017         return ID;
00018     }
00019     void overrideFaceId( int i )
00020     {
00021         ID = i;
00022     }  //USE AT YOUR OWN RISK
00023 
00024   private:
00025     int ID;
00026 };
00027 
00028 //----------------------------------------------------------------------------
00029 typedef std::vector< FaceCellSet > FaceCellSets;
00030 
00031 //----------------------------------------------------------------------------
00032 //class that store the regions that a faces is adjacent too
00033 struct FacesAdjRegions
00034 {
00035     FacesAdjRegions( int f, smoab::EntityHandle r0, smoab::EntityHandle r1 ) : FaceId( f ), Region0( r0 ), Region1( r1 )
00036     {
00037         if( r0 > r1 )
00038         {
00039             std::swap( this->Region0, this->Region1 );
00040         }
00041     }
00042 
00043     FacesAdjRegions( int f ) : FaceId( f ), Region0( -3 ), Region1( -2 ) {}
00044 
00045     bool operator<( const FacesAdjRegions& other ) const
00046     {
00047         return ( this->FaceId < other.FaceId );
00048     }
00049 
00050     smoab::EntityHandle otherId( smoab::EntityHandle other ) const
00051     {
00052         if( other == Region0 )
00053         {
00054             return Region1;
00055         }
00056         return Region0;
00057     }
00058 
00059     int FaceId;
00060     smoab::EntityHandle Region0;
00061     smoab::EntityHandle Region1;
00062 };
00063 //----------------------------------------------------------------------------
00064 smoab::FaceCellSets findFaceSets( smoab::CellSets shells,
00065                                   smoab::CellSets boundaries,
00066                                   std::set< smoab::FacesAdjRegions >& faceMaps )
00067 {
00068     typedef smoab::CellSets::iterator iterator;
00069     typedef smoab::FaceCellSets::iterator faceIterator;
00070     typedef std::set< smoab::FacesAdjRegions >::const_iterator FaceAdjIterator;
00071 
00072     //we need to properly label each unique face in shells
00073     //we do this by intersecting each shell with each other shell
00074     //to find shell on shell contact, and than we intersect each
00075     //resulting shell with the boundary conditions
00076     //the end result of these intersections will be the new modelfaces
00077     int faceId = 1;
00078     smoab::FaceCellSets shellFaces;
00079 
00080     //first intersect each shell with each other shell
00081     std::set< smoab::FacesAdjRegions > shellFaceContacts;
00082     for( iterator i = shells.begin(); i != shells.end(); ++i )
00083     {
00084         //copy the cells so we can add a face that represents
00085         //all the cells of the region that aren't shared with another region
00086         int numCells = i->cells().size();  //size() on range is slow, so cache it
00087         for( iterator j = i + 1; j != shells.end() && numCells > 0; ++j )
00088         {
00089             //intersect i and j to make a new face
00090             smoab::Range intersection = smoab::intersect( i->cells(), j->cells() );
00091             if( !intersection.empty() )
00092             {
00093                 //don't want to increment faceId when the intersection is empty
00094                 smoab::FaceCellSet face( faceId++, i->entity(), intersection );
00095                 shellFaces.push_back( face );
00096                 i->erase( intersection );
00097                 j->erase( intersection );
00098                 numCells -= intersection.size();
00099 
00100                 //add this to the face map
00101                 smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), j->entity() );
00102                 shellFaceContacts.insert( faceInfo );
00103             }
00104         }
00105         //if all the cells for shell i are used, don't add a new
00106         //empty face
00107         if( numCells > 0 )
00108         {
00109             smoab::FaceCellSet face( faceId++, i->entity(), i->cells() );
00110             shellFaces.push_back( face );
00111 
00112             //add this to the face map
00113             smoab::FacesAdjRegions faceInfo( faceId - 1, -1, i->entity() );
00114             shellFaceContacts.insert( faceInfo );
00115         }
00116     }
00117 
00118     //now we have all the faces that match shell on shell contact
00119     //we know process all the new faces to see if they intersect
00120     //with any boundary sets. A boundary set can span multiple
00121     //shells so we want to process it as a second loop
00122 
00123     //store the end before we start adding boundary faces, which
00124     //we don't need to check agianst other boundaries
00125     faceId = 1;  //reset the faced id
00126 
00127     //store in a new face set, expanding the current one causes incorrect results
00128     smoab::FaceCellSets faces;
00129     for( faceIterator i = shellFaces.begin(); i != shellFaces.end(); ++i )
00130     {
00131         //determine from the shell faces if the new face we are creating
00132         //is bounded by two regions or just one
00133         smoab::FacesAdjRegions idToSearchFor( i->faceId() );
00134         FaceAdjIterator adjRegions        = shellFaceContacts.find( idToSearchFor );
00135         smoab::EntityHandle otherRegionId = adjRegions->otherId( i->entity() );
00136 
00137         int numCells = i->cells().size();  //size() on range is slow, so cache it
00138         for( iterator j = boundaries.begin(); j != boundaries.end(); ++j )
00139         {
00140             smoab::Range intersect = smoab::intersect( i->cells(), j->cells() );
00141             if( !intersect.empty() )
00142             {
00143                 //don't want to increment faceId when the intersection is empty
00144                 smoab::FaceCellSet face( faceId++, j->entity(), intersect );
00145                 faces.push_back( face );
00146                 i->erase( intersect );
00147                 numCells -= intersect.size();
00148                 smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), otherRegionId );
00149                 faceMaps.insert( faceInfo );
00150             }
00151         }
00152         if( numCells > 0 )
00153         {
00154             smoab::FaceCellSet face( faceId++, i->entity(), i->cells() );
00155             faces.push_back( face );
00156             smoab::FacesAdjRegions faceInfo( faceId - 1, i->entity(), otherRegionId );
00157             faceMaps.insert( faceInfo );
00158         }
00159     }
00160     return faces;
00161 }
00162 
00163 //----------------------------------------------------------------------------
00164 template < typename T >
00165 std::vector< T > faceIdsPerCell( const smoab::FaceCellSets& faces )
00166 {
00167     typedef smoab::FaceCellSets::const_iterator iterator;
00168     typedef std::vector< smoab::EntityHandle >::const_iterator EntityHandleIterator;
00169     typedef smoab::Range::const_iterator RangeIterator;
00170 
00171     //find all the cells that are in the faceCellSet, and than map
00172     //the proper face id to that relative position, here comes lower_bounds!
00173     std::vector< smoab::EntityHandle > searchableCells;
00174     smoab::Range faceRange = smoab::getAllCells( faces );
00175     smoab::RangeToVector( faceRange, searchableCells );
00176 
00177     //faceIds will be the resulting array
00178     std::vector< T > faceIds( searchableCells.size() );
00179 
00180     //construct the start and end iterators for the lower bounds call
00181 
00182     EntityHandleIterator s_begin = searchableCells.begin();
00183     EntityHandleIterator s_end   = searchableCells.end();
00184 
00185     //search the face cell sets
00186     for( iterator i = faces.begin(); i != faces.end(); ++i )
00187     {
00188         T value                           = static_cast< T >( i->faceId() );
00189         const smoab::Range& entitiesCells = i->cells();
00190         for( RangeIterator j = entitiesCells.begin(); j != entitiesCells.end(); ++j )
00191         {
00192             EntityHandleIterator result = std::lower_bound( s_begin, s_end, *j );
00193             std::size_t newId           = std::distance( s_begin, result );
00194             faceIds[newId]              = value;
00195         }
00196     }
00197     return faceIds;
00198 }
00199 
00200 //----------------------------------------------------------------------------
00201 //given a face adjacency, determine the regions spare tag values
00202 template < typename T >
00203 std::pair< T, T > FaceAdjRegionValues( const smoab::FacesAdjRegions& faceAdj,
00204                                        smoab::Tag* t,
00205                                        const smoab::Interface& interface )
00206 {
00207     std::pair< T, T > returnValue;
00208     const int defaultValue = interface.getDefaultTagVaue< int >( *t );
00209 
00210     /*
00211    *  IF A REGION IS SET TO -1 WE NEED TO PUSH THAT VALUE DOWN
00212    *  AS THE MATERIAL, SINCE THE MOAB DEFAULT TAG VALUE WILL
00213    *  BE CONSIDIERED A REGION, AND WE WANT TO SAY IT BOUNDS THE
00214    *  VOID REGION
00215    */
00216     int tagValue = defaultValue;  //use tagValue to pass in default value
00217     if( faceAdj.Region0 != -1 )
00218     {
00219         tagValue = interface.getTagData( *t, faceAdj.Region0, tagValue );
00220     }
00221     else
00222     {
00223         tagValue = -1;
00224     }
00225 
00226     //set the first region tag value into the pair we are returing
00227     returnValue.first = static_cast< T >( tagValue );
00228 
00229     tagValue = defaultValue;  //use tagValue to pass in default value
00230     tagValue = interface.getTagData( *t, faceAdj.Region1, tagValue );
00231     if( faceAdj.Region1 != -1 )
00232     {
00233         tagValue = interface.getTagData( *t, faceAdj.Region1, tagValue );
00234     }
00235     else
00236     {
00237         tagValue = -1;
00238     }
00239     returnValue.second = static_cast< T >( tagValue );
00240 
00241     return returnValue;
00242 }
00243 
00244 }  // namespace smoab
00245 
00246 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines