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