Branch data Line data Source code
1 : : /*
2 : : * MBVolOp.cpp
3 : : *
4 : : * Created on: Jan 13, 2012
5 : : */
6 : :
7 : : #include "meshkit/MBVolOp.hpp"
8 : : #include "meshkit/MBSplitOp.hpp"
9 : : #include "meshkit/ModelEnt.hpp"
10 : : #include "moab/CartVect.hpp"
11 : :
12 : : //#include "moab/Core.hpp"
13 : : //#include "meshkit/FBiGeom.hpp"
14 : :
15 : :
16 : : namespace MeshKit {
17 : :
18 : :
19 : : //Entity Type initialization for splitting; no mesh output
20 : : moab::EntityType MBVolOp_tps[] = { moab::MBMAXTYPE }; // no mesh, really
21 : 40 : const moab::EntityType* MBVolOp::output_types()
22 : : {
23 : 40 : return MBVolOp_tps;
24 : : }
25 : :
26 : 1 : MBVolOp::MBVolOp(MKCore *mk_core, const MEntVector &me_vec) :
27 [ + - ][ + - ]: 1 : MeshScheme(mk_core, me_vec)
[ + - ]
28 : : {
29 : 1 : _direction[0]=_direction[1]=0.;
30 : 1 : _direction[2]=1.0;
31 : 1 : _pGTT = NULL;
32 : 1 : _rootSet = 0; // it means no set yet, although 0 mean root set in moab
33 : 1 : _fbe = NULL;
34 : 1 : }
35 : :
36 : 0 : MBVolOp::~MBVolOp() {
37 : : // TODO Auto-generated destructor stub
38 [ # # ]: 0 : }
39 : 1 : void MBVolOp::setup_this()
40 : : {
41 : : // construct an FBEngine object, used more like a container, and to trigger the weaving
42 : : // it is not really required; this object is based on gtt object build from top and bottom faces
43 : : // ( which are extracted from model entities of dimension 2)
44 : : // so, involve FBEngine just because we have something we need there
45 : : // collect the top and bottom faces from ment vector
46 : 1 : int nTotSurf = this->mentSelection.size();
47 [ + - ][ + - ]: 1 : std::cout << " total number of faces:" << nTotSurf << "\n";
[ + - ]
48 : :
49 : : // grab all surfaces, and duplicate model using gtt, to get new gsets, that will be
50 : : // continued with volume sets in the end
51 : : // establish the loops from faces
52 [ + - ]: 1 : MEntSelection::iterator mit;
53 : :
54 [ + - ][ + - ]: 1 : moab::Interface * MBI = mk_core()->moab_instance();
55 [ + - ]: 1 : moab::GeomTopoTool gtt(MBI, true);// to retrieve the gsets
56 : : moab::EntityHandle mset;
57 : : // these should all be of dimension 2, faces
58 [ + - ]: 2 : std::vector<moab::EntityHandle> vSurfaces;
59 : :
60 : : moab::ErrorCode rval;
61 : :
62 [ + - ][ + - ]: 7 : for (mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
[ + + ]
63 [ + - ][ + - ]: 6 : mset = (mit->first)->mesh_handle();
64 : : // get the globalId tag
65 [ + - ]: 6 : vSurfaces.push_back(mset);
66 : : }
67 : :
68 [ + - ]: 1 : rval = gtt.duplicate_model(_pGTT, &vSurfaces);
69 [ - + ][ # # ]: 1 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
70 : : // now, _pGTT will have the new sets, which will form the basis for the volume creation
71 : : // new gsets will be added to this _pGTT, and this will (unfortunately) create an OBB tree too,
72 : : // although we do not really need it (actually, 2 obb trees that are not needed!) bummer!
73 : : // the result will be the model root set of the new GTT, stored in result range
74 : : // that corresponds to the first model ent
75 [ + - ]: 1 : ModelEnt * firstMe = (*(mentSelection.begin()) ).first;
76 [ + - ]: 1 : moab::Range & resultRange = mentSelection[firstMe];
77 [ + - ]: 1 : _rootSet = _pGTT->get_root_model_set();
78 [ + - ]: 1 : resultRange.insert(_rootSet);
79 [ + - ][ + - ]: 1 : _fbe = new moab::FBEngine(MBI, _pGTT, false);// smooth=false, not needed here (are you sure?)
80 : : //we will use FBEngine for querying; although ModelEnt would have helped
81 : : // still, the model ent works with geometry adjacency; in this case, there is no geometry
82 : : // we need to work directly with MOAB;
83 : : // not pretty :(
84 : : // we build all this infrastructure and I am not using it
85 : :
86 : :
87 [ + - ]: 1 : establish_mapping();
88 : 2 : return;
89 : : }
90 : 1 : void MBVolOp::establish_mapping()
91 : : {
92 : : // the first half of the surfaces are bottom, next are top
93 : : // establish the connection between top and bottom entities
94 : : // use the direction for vertices, and tangent direction for edges
95 : :
96 : : // these are all vertices from current gtt
97 [ + - ][ + - ]: 1 : moab::Interface * MBI = mk_core()->moab_instance();
98 [ + - ][ + - ]: 1 : moab::Range verticeSets= _pGTT->geoRanges()[0];
99 [ + - ]: 1 : int num_vertices = verticeSets.size();
100 [ + - ]: 2 : std::vector<moab::CartVect> coordVert;
101 [ + - ]: 1 : coordVert.resize(num_vertices);
102 [ + - ]: 2 : std::vector<int> corrV;
103 [ + - ]: 1 : corrV.resize(num_vertices);
104 : 1 : moab::ErrorCode rval = moab::MB_SUCCESS;
105 [ + - ]: 2 : std::map<moab::EntityHandle, int> indexVertex;
106 [ + + ]: 9 : for (int i=0; i<num_vertices; i++)
107 : : {
108 [ + - ][ + - ]: 8 : rval = _fbe->getVtxCoord(verticeSets[i], &(coordVert[i][0]),
109 [ + - ][ + - ]: 16 : &(coordVert[i][1]), &(coordVert[i][2]));
[ + - ][ + - ]
[ + - ][ + - ]
110 [ - + ][ # # ]: 8 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
111 [ + - ]: 8 : corrV[i]=-1;// no correspondence yet
112 [ + - ][ + - ]: 8 : indexVertex[verticeSets[i]] = i;
113 : : }
114 : :
115 : : // decide with an expensive linear search, what vertices correspond to what
116 [ + - ]: 1 : moab::CartVect dirr(_direction);
117 [ + - ]: 1 : dirr.normalize();
118 : 1 : int j=-1;
119 [ + + ]: 9 : for ( j=0; j<num_vertices; j++)
120 : : {
121 [ + - ][ + + ]: 8 : if (corrV[j]!=-1)
122 : 4 : continue; // we already know about this one
123 [ + - ]: 4 : moab::CartVect & v1 = coordVert[j];
124 : 4 : int minIndex = -1;
125 : 4 : double minVal = 1.e38; // HUGE
126 [ + + ]: 36 : for (int k=0; k<num_vertices; k++)
127 : : {
128 [ + + ][ + - ]: 32 : if (j==k || corrV[k]>-1)
[ + + ][ + + ]
129 : 16 : continue;
130 [ + - ][ + - ]: 16 : moab::CartVect product = (v1-coordVert[k])*dirr;
[ + - ]
131 [ + - ]: 16 : double valProd = product.length_squared();
132 [ + + ]: 16 : if (valProd<minVal)
133 : : {
134 : 8 : minVal = valProd;
135 : 16 : minIndex=k;
136 : : }
137 : : }
138 [ + - ][ + - ]: 4 : std::cout<<"j: "<< j <<" val min:" << minVal << " min index:" << minIndex << "\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
139 [ - + ]: 4 : if (minIndex==-1)
140 : : {
141 [ # # ][ # # ]: 0 : std::cout<< "Error Stop it. j: " << j << "\n";
[ # # ]
142 : 0 : continue;
143 : : }
144 : : // make sure that the lower index is on bottom
145 [ + - ][ + - ]: 4 : double dotProd = (coordVert[minIndex]-coordVert[j])%dirr;
[ + - ][ + - ]
146 [ - + ]: 4 : if (dotProd < 0)
147 : : {
148 [ # # ]: 0 : std::cout<<"wrong orientation, bottom vertices should be of lower index\n";
149 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
[ # # ]
150 : : }
151 : :
152 [ + - ]: 4 : corrV[j] = minIndex;
153 [ + - ]: 4 : corrV[minIndex] = j;
154 : : }
155 : : // check that we do not have any left vertices
156 [ + + ]: 9 : for (j=0; j<num_vertices; j++)
157 : : {
158 [ + - ][ - + ]: 8 : if (corrV[j]==-1)
159 : : {
160 [ # # ][ # # ]: 0 : std::cout<<"Error in finding a corresponding vertex to j:"<<j <<"\n";
[ # # ]
161 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
[ # # ]
162 : : }
163 [ + - ][ + - ]: 8 : vertexMap[verticeSets[j]]= verticeSets[corrV[j]];
[ + - ][ + - ]
164 : : }
165 : :
166 : : // so now we have mappings between vertices
167 : : // we need to build mappings between edges and faces.
168 : : // we usually start with bottom and continue to the top
169 [ + - ][ + - ]: 2 : moab::Range edgeSets= _pGTT->geoRanges()[1];
170 [ + - ]: 1 : int numEdges = edgeSets.size();
171 [ + - ]: 2 : std::vector<moab::CartVect> edgeTangents;
172 [ + - ]: 1 : edgeTangents.resize(numEdges*2);// we need 2 tangents per edge
173 : : // for each vertex, compute the tangents at ends, projected on a plan perpendicular to the direction
174 : : // (usually, xy plan...)
175 : : // tangents at both ends!// we know that the u range is from 0 to 1
176 [ + - ]: 2 : std::vector<int> corrE;
177 [ + - ]: 1 : corrE.resize(numEdges);
178 [ + - ]: 2 : std::map<moab::EntityHandle, int> indexEdge;
179 [ + + ]: 13 : for (j=0; j<numEdges; j++)
180 : : {
181 [ + - ]: 12 : corrE[j] = -1;// no correspondence yet
182 [ + - ][ + - ]: 12 : indexEdge[edgeSets[j]] = j;
183 : : // careful about the FBEngine, it is not smooth
184 [ + - ]: 12 : std::vector<moab::EntityHandle> meshEdges;
185 [ + - ][ + - ]: 12 : rval = MBI->get_entities_by_type(edgeSets[j], moab::MBEDGE, meshEdges);
186 [ - + ][ # # ]: 12 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
187 : : // get first and last edges, connectivity and compute tangent
188 : : // it will be used to map edges
189 : : // first edge, last edge
190 [ + - ]: 12 : moab::EntityHandle firstMeshEdge = meshEdges[0];
191 [ + - ]: 12 : moab::EntityHandle lastMeshEdge = meshEdges[meshEdges.size()-1];
192 : 12 : const moab::EntityHandle * conn=NULL;
193 : : int nnodes;
194 [ + - ]: 12 : rval = MBI->get_connectivity(firstMeshEdge, conn, nnodes);
195 [ - + ][ # # ]: 12 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
196 [ + - ][ + + ]: 36 : moab::CartVect posi[2];
197 [ + - ][ + - ]: 12 : rval = MBI->get_coords(conn, 2, &(posi[0][0]));
198 [ - + ][ # # ]: 12 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
199 [ + - ]: 12 : posi[0] = posi[1]-posi[0]; //
200 [ + - ]: 12 : posi[0] = posi[0]*dirr;
201 [ + - ][ + - ]: 12 : edgeTangents[2*j] = posi[0]*dirr;// this is in the plane
202 [ + - ][ + - ]: 12 : edgeTangents[2*j].normalize();// it should be non null, but accidents happen :)
203 [ + - ]: 12 : rval = MBI->get_connectivity(lastMeshEdge, conn, nnodes);
204 [ - + ][ # # ]: 12 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
205 [ + - ][ + - ]: 12 : rval = MBI->get_coords(conn, 2, &(posi[0][0]));
206 [ - + ][ # # ]: 12 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
207 [ + - ]: 12 : posi[0] = posi[1]-posi[0]; //
208 [ + - ]: 12 : posi[0] = posi[0]*dirr;
209 [ + - ][ + - ]: 12 : edgeTangents[2*j+1] = posi[0]*dirr;// this is in the plane
210 [ + - ][ + - ]: 12 : edgeTangents[2*j+1].normalize();
211 : 12 : }
212 : : // now try to match edges based on their vertices matching, and start and end tangents matching
213 [ + + ]: 13 : for (j = 0; j<numEdges; j++)
214 : : {
215 [ + - ][ + + ]: 12 : if (corrE[j]>=0)
216 : 6 : continue; // we already have a correspondent edge for it
217 [ + - ]: 6 : moab::Range adjVertices;
218 [ + - ][ + - ]: 6 : rval = _fbe-> getEntAdj(edgeSets[j], /*vertex type*/0, adjVertices);
219 [ - + ][ # # ]: 6 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
220 [ + - ][ + - ]: 6 : if (adjVertices.size()==2)
221 : : {
222 : 6 : int sense=0;
223 [ + - ][ + - ]: 6 : rval = _fbe->getEgVtxSense(edgeSets[j], adjVertices[0], adjVertices[1], sense);
[ + - ][ + - ]
224 [ - + ][ # # ]: 6 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
225 : : // find the edges that are adjacent to map(v0) and map(v1)
226 [ + - ][ + - ]: 6 : moab::EntityHandle mapV0 = vertexMap [adjVertices[0]];
227 [ + - ][ + - ]: 6 : moab::EntityHandle mapV1 = vertexMap [adjVertices[1]];
228 : : // get all edges adjacent to both vertices
229 [ + - ][ + - ]: 12 : moab::Range adjEdges0, adjEdges1;
230 [ + - ]: 6 : rval = _fbe-> getEntAdj(mapV0, /*edge type*/1, adjEdges0);
231 [ - + ][ # # ]: 6 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
232 [ + - ]: 6 : rval = _fbe-> getEntAdj(mapV1, /*edge type*/1, adjEdges1);
233 [ + - ][ + - ]: 6 : adjEdges0 = intersect(adjEdges0, adjEdges1);
234 : : // the mapped edge should be in this range
235 [ + - ]: 6 : moab::CartVect & t0=edgeTangents[2*j];
236 [ + - ]: 6 : moab::CartVect & t1=edgeTangents[2*j+1];
237 [ + - ][ + - ]: 16 : for (moab::Range::iterator rit= adjEdges0.begin(); rit!=adjEdges0.end(); rit++)
[ + - ][ + - ]
[ + + ]
238 : : {
239 [ + - ]: 10 : moab::EntityHandle candidateMapEdge = *rit;
240 : 10 : int sense1=0;
241 [ + - ]: 10 : int indexCandEdge = indexEdge[candidateMapEdge];
242 [ - + ]: 10 : if (indexCandEdge == j)
243 : : // error
244 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);
[ # # ]
245 [ + - ]: 10 : moab::CartVect & tm0=edgeTangents[2*indexCandEdge];
246 [ + - ]: 10 : moab::CartVect & tm1=edgeTangents[2*indexCandEdge+1];
247 [ + - ]: 10 : rval = _fbe->getEgVtxSense(candidateMapEdge, mapV0, mapV1, sense1);
248 [ - + ][ # # ]: 10 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
249 : :
250 [ + + ]: 10 : if (sense*sense1>0)// same sense
251 : : {
252 [ + - ][ + + ]: 8 : if ( (t0%tm0 >= 0.99 ) && (t1%tm1>=0.99))// same sense, almost..
[ + - ][ + - ]
[ + + ]
253 : : {
254 [ + - ]: 6 : corrE[j]=indexCandEdge;
255 [ + - ]: 6 : corrE[indexCandEdge] = j;
256 : : }
257 : : }
258 : : else
259 : : {
260 : : // different senses, check the opposite tangents
261 [ + - ][ - + ]: 2 : if ( (t0%tm0 <=-0.99) && (t1%tm1<=-0.99) )
[ # # ][ # # ]
[ - + ]
262 : : {
263 : : // how do we store the opposite senses?
264 [ # # ]: 0 : corrE[j]=indexCandEdge;
265 [ # # ]: 0 : corrE[indexCandEdge] = j;
266 : : }
267 : : }
268 : : }
269 [ + - ][ - + ]: 6 : if (corrE[j]==-1)
270 [ # # ][ # # ]: 6 : MBERRCHK(moab::MB_FAILURE, MBI);// can't find a corresponding edge
[ # # ]
271 : :
272 : : }
273 : 6 : }
274 : : // so we have edge matching for every edge
275 [ + + ]: 13 : for (j=0; j<numEdges; j++)
276 : : {
277 [ + - ][ - + ]: 12 : if (corrE[j]==-1)
278 : : {
279 [ # # ][ # # ]: 0 : std::cout<<"Error in finding a corresponding edge: "<<j <<"\n";
[ # # ]
280 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
[ # # ]
281 : : }
282 [ + - ][ + - ]: 12 : std::cout<< "edge j=" << j << " mapped to edge " << corrE[j] << "\n";
[ + - ][ + - ]
[ + - ][ + - ]
283 [ + - ][ + - ]: 12 : edgeMap[edgeSets[j]]= edgeSets[corrE[j]];
[ + - ][ + - ]
284 : : }
285 : : // now, we still need to separate top and bottom faces, somehow
286 : : // we assume faces were stenciled correctly with the same polylines
287 : : // start from bottom vertices (first half of the vertices is on bottom)
288 [ + - ]: 2 : moab::Range bottomFaces;
289 [ + + ]: 5 : for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
290 : : {
291 [ + - ]: 4 : moab::Range faces;
292 [ + - ][ + - ]: 4 : rval = _fbe-> getEntAdj(verticeSets[j], /*face type*/2, faces);
293 [ - + ][ # # ]: 4 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
294 [ + - ]: 4 : bottomFaces.merge(faces);
295 : 4 : }
296 : : // now find the mapped face for each of these, based on mapped edges
297 [ + - ][ + + ]: 4 : for (j=0; j<(int)bottomFaces.size(); j++)
298 : : {
299 [ + - ]: 3 : moab::Range edges;
300 [ + - ][ + - ]: 3 : rval = _fbe-> getEntAdj(bottomFaces[j], /*edge type*/1, edges);
301 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
302 [ + - ]: 6 : moab::Range mappedEdges;
303 : : // get now the face connected to the mapped edges
304 : 3 : int i=0;
305 : : // find the face with all those edges among the other faces
306 : :
307 [ + - ]: 6 : moab::Range mapFaces;
308 [ + - ][ + - ]: 3 : rval = _fbe-> getEntAdj(edgeMap[edges[0]], /*edge type*/2, mapFaces);
[ + - ]
309 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
310 [ + - ][ + + ]: 8 : for (i=1; i<(int)edges.size(); i++)
311 : : {
312 [ + - ]: 5 : moab::Range faces;
313 [ + - ][ + - ]: 5 : rval = _fbe-> getEntAdj(edgeMap[edges[i]], /*edge type*/2, faces);
[ + - ]
314 [ + - ][ + - ]: 5 : mapFaces = intersect(mapFaces, faces);
315 : 5 : }
316 [ + - ][ - + ]: 3 : if (mapFaces.size()!=1)
317 : : {
318 [ # # ][ # # ]: 0 : std::cout<<"Can't find unique mapped face to face index " << j << "\n";
[ # # ]
319 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
[ # # ]
320 : : }
321 [ + - ][ + - ]: 3 : faceMap[bottomFaces[j]] = mapFaces[0];
[ + - ]
322 [ + - ][ + - ]: 3 : faceMap[mapFaces[0]]=bottomFaces[j];
[ + - ]
323 [ + - ][ + - ]: 3 : std::cout << "face j:" << j << " set:" << MBI->id_from_handle(bottomFaces[j]) << " mapped to "
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
324 [ + - ][ + - ]: 6 : << MBI->id_from_handle(mapFaces[0]) << "\n";
[ + - ][ + - ]
325 : 4 : }
326 : :
327 : 1 : }
328 : :
329 : 1 : void MBVolOp::execute_this()
330 : : {
331 : : // now, we have the maps between top and bottom faces established, also for edges and vertices
332 : : // first build edges between corresponding vertices, then faces between edges, then
333 : : // finally, volumes
334 : : // we start from bottom towards top
335 : 1 : moab::ErrorCode rval = moab::MB_SUCCESS;
336 [ + - ][ + - ]: 1 : moab::Interface * MBI = mk_core()->moab_instance();
337 : :
338 [ + - ][ + - ]: 1 : moab::Range verticeSets=_pGTT->geoRanges()[0];
339 [ + - ]: 1 : int num_vertices = verticeSets.size();
340 [ + - ]: 2 : moab::Range bottomEdges;
341 : 1 : int j=0;
342 [ + + ]: 5 : for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
343 : : {
344 [ + - ]: 4 : moab::Range edges;
345 [ + - ][ + - ]: 4 : rval = _fbe-> getEntAdj(verticeSets[j], /*edge type*/1, edges);
346 [ - + ][ # # ]: 4 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
347 [ + - ]: 4 : bottomEdges.merge(edges);
348 : 4 : }
349 : : // we need to decide bottom faces before we create new faces, edges, etc
350 [ + - ]: 2 : moab::Range bottomFaces;
351 [ + + ]: 5 : for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
352 : : {
353 [ + - ]: 4 : moab::Range faces;
354 [ + - ][ + - ]: 4 : rval = _fbe-> getEntAdj(verticeSets[j], /*face type*/2, faces);
355 [ - + ][ # # ]: 4 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
356 [ + - ]: 4 : bottomFaces.merge(faces);
357 : 4 : }
358 : :
359 : : // we have verified that the first half are from bottom surfaces
360 : : // for each bottom face we will create a volume;
361 : : // start with the edges, create weaving faces between edges (no new vertices!)
362 : : // now, for each bottom edge, create a weaving face
363 [ + - ]: 2 : std::vector<moab::EntityHandle> newFaces;
364 [ + - ][ + - ]: 1 : newFaces.resize(bottomEdges.size());
365 [ + - ][ + + ]: 7 : for (j=0; j<(int)bottomEdges.size(); j++)
366 : : {
367 [ + - ][ + - ]: 6 : rval = _fbe->weave_lateral_face_from_edges(bottomEdges[j], edgeMap[bottomEdges[j]], _direction, newFaces[j]);
[ + - ][ + - ]
[ + - ]
368 [ - + ][ # # ]: 6 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
369 : : }
370 : : // now, to create volumes, we need to loop over bottom faces and add volumes one by one, and the
371 : : // orientation of faces in them
372 : :
373 : 1 : int volumeMatId = 0;// just give a mat id, for debugging, mostly
374 : : moab::Tag matTag;
375 [ + - ]: 1 : rval = MBI->tag_get_handle("MATERIAL_SET", 1, moab::MB_TYPE_INTEGER, matTag);
376 [ - + ][ # # ]: 1 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
377 : :
378 [ + - ][ + + ]: 4 : for (j = 0; j<(int)bottomFaces.size(); j++)
379 : : {
380 : : // create volume here , finally!!!
381 : : moab::EntityHandle volume;
382 [ + - ]: 3 : rval = MBI->create_meshset(moab::MESHSET_SET, volume);
383 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
384 : 3 : volumeMatId++;
385 [ + - ]: 3 : rval= MBI->tag_set_data(matTag, &volume, 1, &volumeMatId);
386 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
387 [ + - ]: 3 : moab::EntityHandle botFace = bottomFaces[j];
388 [ + - ]: 3 : moab::EntityHandle topFace = faceMap[botFace];
389 : : // start copy
390 : : // get the edges of bot face
391 [ + - ]: 3 : rval= MBI->add_parent_child(volume, botFace);
392 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
393 : :
394 [ + - ]: 3 : rval= MBI->add_parent_child(volume, topFace);
395 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
396 : :
397 [ + - ]: 3 : rval = _pGTT->add_geo_set(volume, 3);
398 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
399 : :
400 : : // set senses
401 : : // bottom face is negatively oriented, its normal is toward interior of the volume
402 [ + - ]: 3 : rval = _pGTT->set_sense(botFace, volume, -1);
403 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
404 : :
405 : : // the top face is positively oriented
406 [ + - ]: 3 : rval = _pGTT->set_sense(topFace, volume, 1);
407 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
408 : :
409 : : // the children should be in the same direction
410 : : // get the side edges of each face, and form lateral faces, along direction
411 [ + - ]: 3 : std::vector<moab::EntityHandle> edges1;
412 : :
413 [ + - ]: 3 : rval = MBI->get_child_meshsets(botFace, edges1); // no hops
414 [ - + ][ # # ]: 3 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
415 : :
416 [ + + ]: 11 : for (unsigned int i = 0; i < edges1.size(); ++i)
417 : : {
418 : : // the orientation of edge in face will give the sense of face in volume
419 [ + - ][ + - ]: 8 : int indexB = bottomEdges.index(edges1[i]);
420 [ - + ]: 8 : if (indexB<=-1)
421 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, MBI);
[ # # ]
422 : 8 : int sense = 1;
423 [ + - ][ + - ]: 8 : rval = _pGTT->get_sense(edges1[i], botFace, sense);
424 [ - + ][ # # ]: 8 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
425 [ + - ][ + - ]: 8 : rval=MBI->add_parent_child(volume, newFaces[indexB]);
426 [ - + ][ # # ]: 8 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
427 : :
428 : : // set sense as the one decided from the bottom edge sense within the bottom face
429 [ + - ][ + - ]: 8 : rval = _pGTT->set_sense(newFaces[indexB], volume, sense);
430 [ - + ][ # # ]: 8 : MBERRCHK(rval, MBI);
[ # # ][ # # ]
431 : : }
432 : : // end copy
433 : 3 : }
434 : :
435 : :
436 [ + - ]: 1 : delete _fbe;
437 [ + - ]: 2 : delete _pGTT; // when we are done, remove the _pGTT;
438 : : // at this point, the result will be the model root set of the first ment of the model
439 : : // ( (*(mentSelection.begin()) ).second )
440 : 1 : }
441 : :
442 [ + - ][ + - ]: 156 : }
|
