MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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