Branch data Line data Source code
1 : : #include "meshkit/TFIMapping.hpp"
2 : : #include "meshkit/MKCore.hpp"
3 : : #include "meshkit/EdgeMesher.hpp"
4 : : #include "meshkit/ModelEnt.hpp"
5 : : #include "meshkit/SizingFunction.hpp"
6 : : #include "moab/ReadUtilIface.hpp"
7 : : #include "EquipotentialSmooth.hpp"
8 : : #ifdef HAVE_MESQUITE
9 : : #include "meshkit/MeshImprove.hpp"
10 : : #endif
11 : :
12 : : #include <vector>
13 : : #include <iostream>
14 : : #include <math.h>
15 : : #include <map>
16 : : #include <algorithm>
17 : :
18 : :
19 : : namespace MeshKit {
20 : :
21 : : //---------------------------------------------------------------------------//
22 : : //Entity Type initialization for TFIMapping meshing
23 : : moab::EntityType TFIMapping_tps[] = { moab::MBVERTEX, moab::MBEDGE, moab::MBQUAD, moab::MBMAXTYPE };
24 : 42 : const moab::EntityType* TFIMapping::output_types()
25 : : {
26 : 42 : return TFIMapping_tps;
27 : : }
28 : :
29 : : //---------------------------------------------------------------------------//
30 : : // construction function for TFIMapping class
31 : 7 : TFIMapping::TFIMapping(MKCore *mk_core, const MEntVector &me_vec) :
32 : 7 : MeshScheme(mk_core, me_vec)
33 : : {
34 : 7 : _shapeImprove=false;
35 : 7 : }
36 : :
37 : : //---------------------------------------------------------------------------//
38 : : // deconstruction function for TFIMapping class
39 : 21 : TFIMapping::~TFIMapping()
40 : : {
41 : :
42 [ - + ]: 14 : }
43 : :
44 : : //---------------------------------------------------------------------------//
45 : : // setup function:
46 : 7 : void TFIMapping::setup_this()
47 : : {
48 : :
49 : : /* the only things we need to make sure :
50 : : 1) there are 4 edges, exactly: this is removed
51 : : 2) the opposite edges have the same meshcount
52 : : - if some of the edges are meshed, we need to mesh the opposite edge with
53 : : correct meshcount
54 : : - if 2 opposite edges are meshed, verify the mesh count
55 : : */
56 : : // get iGeom instance from the first ment selection
57 [ - + ]: 7 : if (mentSelection.empty())
58 : 7 : return;
59 : :
60 : : //loop over the surfaces
61 [ + - ][ + - ]: 32 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
62 : : {
63 [ + - ]: 25 : ModelEnt *me = mit -> first;
64 [ + - ]: 25 : int dimME = me->dimension();
65 [ - + ]: 25 : if (dimME != 2)
66 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "bad input for TFI Mapping, we can only mesh surfaces");
67 : : //first check whether the surface is meshed or not
68 [ + - ][ - + ]: 25 : if (me->get_meshed_state() >= COMPLETE_MESH)
69 : 0 : continue;
70 [ + - ]: 25 : int size_index = me->sizing_function_index();
71 : 25 : int mesh_count_surface = -1;
72 [ + - ]: 25 : if (size_index >= 0)
73 [ + - ][ + - ]: 25 : mesh_count_surface = me->mk_core()->sizing_function(size_index)->intervals();
[ + - ]
74 : :
75 : : // get the boundary loop of edges
76 [ + - ]: 25 : MEntVector boundEdges;
77 [ + - ][ + - ]: 50 : std::vector<int> senses, group_sizes;
78 [ + - ]: 25 : me->ModelEnt::boundary(1, boundEdges, &senses, &group_sizes);
79 : : //remove this constraint in case of the side-cylinder surface
80 : : //if (boundEdges.size() != 4)
81 : : // ECERRCHK(MK_FAILURE, "bad input for TFI Mapping, we can only mesh surfaces with 4 edges");
82 : :
83 [ - + ]: 25 : if (boundEdges.size()<4){
84 [ # # ][ # # ]: 0 : ModelEnt *oppEdges[2] = {boundEdges[0], boundEdges[1]};
85 [ # # ]: 0 : MEntVector edgesToMesh;
86 : :
87 : 0 : int mesh_count = mesh_count_surface;
88 : 0 : bool force = false;
89 [ # # ]: 0 : for (int j = 0; j < 2; j++){
90 [ # # ][ # # ]: 0 : if (oppEdges[j]->get_meshed_state() >= COMPLETE_MESH){
91 [ # # ]: 0 : std::vector<moab::EntityHandle> medges;
92 [ # # ]: 0 : oppEdges[j]->get_mesh(1, medges, true);
93 : 0 : mesh_count = (int) medges.size();
94 : 0 : force = true;
95 : : }
96 : : else
97 : : {
98 [ # # ]: 0 : int indexS = oppEdges[j]->sizing_function_index();
99 [ # # ]: 0 : if (indexS >= 0){
100 [ # # ][ # # ]: 0 : SizingFunction * sfe = mk_core()->sizing_function(indexS);
101 [ # # ]: 0 : if (!force)
102 : : {
103 [ # # ][ # # ]: 0 : if (sfe->intervals() > 0)
104 [ # # ]: 0 : mesh_count = sfe->intervals();
105 [ # # ][ # # ]: 0 : else if (sfe->size() > 0)
106 [ # # ][ # # ]: 0 : mesh_count = oppEdges[j]->measure() / sfe->size();
107 [ # # ][ # # ]: 0 : if (mesh_count % 2 && oppEdges[j]->constrain_even())
[ # # ][ # # ]
108 : 0 : ++mesh_count;
109 : : }
110 : : }
111 [ # # ]: 0 : edgesToMesh.push_back(oppEdges[j]);
112 : : }
113 : : }
114 [ # # ]: 0 : if (edgesToMesh.size() > 0)
115 : : {
116 [ # # ][ # # ]: 0 : EdgeMesher * em = (EdgeMesher*) me->mk_core()->construct_meshop("EdgeMesher", edgesToMesh);
[ # # ]
117 [ # # ]: 0 : if (mesh_count < 0)
118 : : {
119 [ # # ]: 0 : std::cout << "mesh count not set properly on opposite edges, set it to 10\n";
120 : 0 : mesh_count = 10; // 4 is a nice number, used in the default edge mesher;
121 : : // but I like 10 more
122 : : }
123 : :
124 [ # # ]: 0 : for (unsigned int j = 0; j < edgesToMesh.size(); j++)
125 : : {
126 : 0 : int edgeMeshCount = 0;
127 : : int edgeSfIndex =
128 [ # # ][ # # ]: 0 : edgesToMesh[j]->sizing_function_index();
129 [ # # ]: 0 : if (edgeSfIndex >= 0)
130 : : {
131 : : SizingFunction* edgeSf =
132 [ # # ][ # # ]: 0 : mk_core()->sizing_function(edgeSfIndex);
133 [ # # ]: 0 : edgeMeshCount = edgeSf->intervals();
134 : : }
135 [ # # ]: 0 : if (mesh_count != edgeMeshCount)
136 : : {
137 [ # # ][ # # ]: 0 : edgesToMesh[j]->mesh_intervals(mesh_count);
138 : : }
139 [ # # ]: 0 : if (force)
140 : : {
141 : : // the opposite edge is already meshed, so the number
142 : : // of intervals is a hard constraint for this edge
143 [ # # ][ # # ]: 0 : edgesToMesh[j]->interval_firmness(HARD);
144 : : }
145 [ # # ][ # # ]: 0 : edgesToMesh[j]->add_meshop(em);
146 : : }
147 [ # # ]: 0 : mk_core()->insert_node(em, (GraphNode*)this,
148 [ # # ][ # # ]: 0 : mk_core()->root_node());
[ # # ]
149 : 0 : }
150 : : }
151 : : else{
152 : :
153 : : // mesh edge 0 and 2 together, and 1 and 3 together (same mesh count)
154 : : // look at all settings, to decide proper mesh count
155 [ + + ]: 75 : for (int k = 0; k <= 1; k++)
156 : : {
157 : : // treat first edges 0 and 2, then 1 and 3
158 [ + - ][ + - ]: 50 : ModelEnt * oppEdges[2] = { boundEdges[k], boundEdges[k + 2] };
159 [ + - ]: 50 : MEntVector edgesToMesh;// edges that are not meshed yet
160 : : //if one of them is meshed and the other not, use the same mesh count
161 : :
162 : : // take the maximum of the proposed mesh counts, either from sizing function, or mesh intervals
163 : 50 : int mesh_count = mesh_count_surface; // could be -1, still
164 : 50 : bool force = false;
165 [ + + ]: 150 : for (int j = 0; j < 2; j++)
166 : : {
167 [ + - ][ + + ]: 100 : if (oppEdges[j]->get_meshed_state() >= COMPLETE_MESH)
168 : : {
169 : : // in this case, force the other edge to have the same mesh count, do not take it from surface
170 [ + - ]: 12 : std::vector<moab::EntityHandle> medges;
171 [ + - ]: 12 : oppEdges[j]->get_mesh(1, medges, true);
172 : 12 : mesh_count = (int) medges.size();
173 : 12 : force = true;
174 : : }
175 : : else
176 : : {
177 [ + - ]: 88 : int indexS = oppEdges[j]->sizing_function_index();
178 [ + - ]: 88 : if (indexS >= 0)
179 : : {
180 [ + - ][ + - ]: 88 : SizingFunction * sfe = mk_core()->sizing_function(indexS);
181 [ + + ]: 88 : if (!force)
182 : : {
183 : : // if a sizing function was set on an edge, use
184 : : // that rather than a mesh count from the surface
185 [ + - ][ + + ]: 82 : if (sfe->intervals() > 0)
186 [ + - ]: 18 : mesh_count = sfe->intervals();
187 [ + - ][ + - ]: 64 : else if (sfe->size() > 0)
188 [ + - ][ + - ]: 64 : mesh_count = oppEdges[j]->measure() / sfe->size();
189 [ + + ][ + - ]: 82 : if (mesh_count % 2 && oppEdges[j]->constrain_even())
[ + + ][ + + ]
190 : 88 : ++mesh_count;
191 : : }
192 : : }
193 : : // push it to the list if it is not setup to another mesh op (edge mesher) already
194 : : //if (oppEdges[j]->is_meshops_list_empty())// it will create an EdgeMesher later
195 [ + + ][ + + ]: 172 : if ((j == 0 || (oppEdges[j] != oppEdges[0])) &&
[ + + ][ + + ]
196 [ + - ]: 84 : oppEdges[j]->is_meshops_list_empty())
197 : : {
198 [ + - ]: 36 : edgesToMesh.push_back(oppEdges[j]);
199 : : }
200 : : }
201 : : }
202 : : // decide on a mesh count now, if edgesToMesh.size()>0
203 [ + + ]: 50 : if (edgesToMesh.size() > 0)
204 : : {
205 : :
206 [ + - ][ + - ]: 30 : EdgeMesher * em = (EdgeMesher*) me->mk_core()->construct_meshop("EdgeMesher", edgesToMesh);
[ + - ]
207 [ - + ]: 30 : if (mesh_count < 0)
208 : : {
209 [ # # ]: 0 : std::cout << "mesh count not set properly on opposite edges, set it to 10\n";
210 : 0 : mesh_count = 10; // 4 is a nice number, used in the default edge mesher;
211 : : // but I like 10 more
212 : : }
213 : :
214 [ + + ]: 66 : for (unsigned int j = 0; j < edgesToMesh.size(); j++)
215 : : {
216 : 36 : int edgeMeshCount = 0;
217 : : int edgeSfIndex =
218 [ + - ][ + - ]: 36 : edgesToMesh[j]->sizing_function_index();
219 [ + - ]: 36 : if (edgeSfIndex >= 0)
220 : : {
221 : : SizingFunction* edgeSf =
222 [ + - ][ + - ]: 36 : mk_core()->sizing_function(edgeSfIndex);
223 [ + - ]: 36 : edgeMeshCount = edgeSf->intervals();
224 : : }
225 [ + + ]: 36 : if (mesh_count != edgeMeshCount)
226 : : {
227 [ + - ][ + - ]: 28 : edgesToMesh[j]->mesh_intervals(mesh_count);
228 : : }
229 [ + + ]: 36 : if (force)
230 : : {
231 : : // the opposite edge is already meshed, so the number
232 : : // of intervals is a hard constraint for this edge
233 [ + - ][ + - ]: 8 : edgesToMesh[j]->interval_firmness(HARD);
234 : : }
235 [ + - ][ + - ]: 36 : edgesToMesh[j]->add_meshop(em);
236 : : }
237 [ + - ]: 30 : mk_core()->insert_node(em, (GraphNode*)this,
238 [ + - ][ + - ]: 60 : mk_core()->root_node());
[ + - ]
239 : : }
240 : 50 : } // end loop over pair of opposite edges
241 : : }
242 : 25 : }// end loop over surfaces
243 : :
244 : 7 : ensure_facet_dependencies(false);
245 : :
246 : 7 : mk_core()->print_graph("AfterTFISetup.eps");
247 : : }
248 : :
249 : : //---------------------------------------------------------------------------//
250 : : // execute function: generate the all-quad mesh through the TFI mapping
251 : 7 : void TFIMapping::execute_this()
252 : : {
253 : :
254 : : //loop over the surfaces
255 [ + - ][ + - ]: 32 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
256 : : {
257 [ + - ]: 25 : ModelEnt *me = mit -> first;
258 : : //first check whether the surface is meshed or not
259 [ + - ][ + + ]: 25 : if (me->get_meshed_state() >= COMPLETE_MESH){
260 : :
261 : : #ifdef HAVE_MESQUITE
262 : : iBase_EntitySetHandle entityset;
263 [ + - ][ + - ]: 6 : iRel::Error r_err = mk_core()->irel_pair(me->iRelPairIndex())->getEntSetRelation(me->geom_handle(), 0, entityset);
[ + - ][ + - ]
[ + - ]
264 [ + - ]: 6 : IBERRCHK(r_err, "Trouble get the entityset w.r.t a surface!");
265 [ + - ][ + - ]: 6 : MeshImprove shapesmooth(mk_core(), false, false, true, false, mk_core()->igeom_instance(me->iGeomIndex()));
[ + - ][ + - ]
[ + - ]
266 [ + - ][ + - ]: 6 : shapesmooth.SurfMeshImprove(me->geom_handle(), entityset, iBase_FACE);
267 : : #endif
268 : :
269 : 6 : continue;
270 : : }
271 : :
272 [ + - ]: 19 : MEntVector boundEdges;
273 [ + - ][ + - ]: 38 : std::vector<int> senses, group_sizes;
274 [ + - ]: 19 : me->ModelEnt::boundary(1, boundEdges, &senses, &group_sizes);
275 [ + - ]: 38 : set<ModelEnt*> distinctBoundEdges;
276 [ + - ]: 19 : distinctBoundEdges.insert(boundEdges.begin(), boundEdges.end());
277 [ + + ]: 19 : if (distinctBoundEdges.size() == 4)
278 [ + - ]: 17 : SurfMapping(me);
279 : : else
280 [ + - ]: 2 : cylinderSurfMapping(me);
281 : :
282 : : //ok, we are done, commit to ME
283 [ + - ][ + - ]: 19 : me->commit_mesh(mit->second, COMPLETE_MESH);
284 : 19 : }
285 : :
286 : 7 : }
287 : :
288 : 2 : int TFIMapping::cylinderSurfMapping(ModelEnt *ent)
289 : : {
290 [ + - ]: 2 : int irelPairIndex = ent->iRelPairIndex();
291 : :
292 : : // determine whether there is an edge along the linking surface
293 [ + - ]: 2 : MEntVector allBoundEdges;
294 [ + - ][ + - ]: 4 : std::vector<int> allSenses, group_sizes;
295 [ + - ]: 2 : ent->ModelEnt::boundary(1, allBoundEdges, &allSenses, &group_sizes);
296 [ + - ]: 4 : map<ModelEnt*, int> boundEdgeCount;
297 [ + + ]: 10 : for (unsigned int i = 0; i < allBoundEdges.size(); ++i)
298 : : {
299 [ + - ][ + - ]: 8 : boundEdgeCount[allBoundEdges[i]] = 0;
300 : : }
301 [ + + ]: 10 : for (unsigned int i = 0; i < allBoundEdges.size(); ++i)
302 : : {
303 [ + - ][ + - ]: 8 : ++boundEdgeCount[allBoundEdges[i]];
304 : : }
305 [ + - ]: 4 : MEntVector boundEdges;
306 [ + - ]: 4 : std::vector<int> boundEdgeSenses;
307 : 2 : ModelEnt* linkSurfEdge = NULL;
308 [ + + ]: 10 : for (unsigned int i = 0; i < allBoundEdges.size(); ++i)
309 : : {
310 [ + - ][ + - ]: 8 : if (boundEdgeCount[allBoundEdges[i]] == 1)
[ + + ]
311 : : {
312 [ + - ][ + - ]: 4 : boundEdges.push_back(allBoundEdges[i]);
313 [ + - ][ + - ]: 4 : boundEdgeSenses.push_back(allSenses[i]);
314 : : }
315 : : else
316 : : {
317 [ + - ]: 4 : linkSurfEdge = allBoundEdges[i];
318 : : }
319 : : }
320 : :
321 [ - + ]: 2 : if (boundEdges.size() != 2)
322 : : {
323 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Cylinder TFIMapping does not have exactly two distinct bounding edges");
324 : : }
325 : :
326 [ + - ]: 2 : std::cout << "TFIMapping on cylinder\n";
327 : :
328 [ + - ]: 4 : std::vector<moab::EntityHandle> nList;
329 [ + - ][ + - ]: 4 : std::vector<iBase_EntityHandle> List_i, List_ii;
330 : :
331 : : // get nodes of bounding edge 0 into List_i
332 [ + - ][ + - ]: 2 : boundEdges[0]->get_mesh(0, nList, true);
333 : 2 : unsigned int ix = 0;
334 : 2 : int size_i = (int)nList.size()-1;
335 [ + + ]: 36 : for (ix = 0; ix < nList.size(); ix++)
336 : : {
337 [ + - ][ + - ]: 34 : List_i.push_back((iBase_EntityHandle) nList[ix]);
338 : : }
339 : : // we reverse the first boundary edge if it is "negative" in the loop
340 [ + - ][ + - ]: 2 : if (boundEdgeSenses[0] == -1)
341 : : {
342 [ + - ]: 2 : std::reverse(List_i.begin(), List_i.end());
343 : : }
344 : 2 : nList.clear();
345 : :
346 : : // get nodes of bounding edge 1 into List_ii
347 [ + - ][ + - ]: 2 : boundEdges[1]->get_mesh(0, nList, true);
348 : 2 : int size_ii = nList.size() - 1;
349 [ + + ]: 36 : for (ix = 0; ix < nList.size(); ix++)
350 : : {
351 [ + - ][ + - ]: 34 : List_ii.push_back((iBase_EntityHandle) nList[ix]);
352 : : }
353 : : // we reverse the second boundary edge if it is "positive" in the loop
354 [ + - ][ + - ]: 2 : if (boundEdgeSenses[1] == 1)
355 : : {
356 [ + - ]: 2 : std::reverse(List_ii.begin(), List_ii.end());
357 : : }
358 : :
359 [ - + ]: 2 : if (size_i != size_ii)
360 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Opposite edges have different mesh count, abort");
361 : :
362 : : // get nodes that are on the linking surface edge, if any,
363 : : // and identify where they start on the source and target
364 : 2 : int linkingEdgeNodeI = -1;
365 : 2 : int linkingEdgeNodeII = -1;
366 : 2 : int offset = 0;
367 [ + - ]: 4 : std::vector<iBase_EntityHandle> linkEdgeNodeList;
368 [ + - ]: 2 : if (linkSurfEdge != NULL)
369 : : {
370 [ + - ]: 2 : std::vector<moab::EntityHandle> lseNodes;
371 [ + - ]: 2 : linkSurfEdge->get_mesh(0, lseNodes, true);
372 [ + + ]: 8 : for (ix = 0; ix < lseNodes.size(); ++ix)
373 : : {
374 [ + - ][ + - ]: 6 : linkEdgeNodeList.push_back((iBase_EntityHandle)lseNodes[ix]);
375 : : }
376 [ + - ]: 2 : iBase_EntityHandle nodeOnI = linkEdgeNodeList[0];
377 [ + - ]: 2 : iBase_EntityHandle nodeOnII = linkEdgeNodeList[linkEdgeNodeList.size() - 1];
378 [ + - ]: 2 : for (ix = 0; ix < List_i.size(); ix++)
379 : : {
380 [ + - ][ - + ]: 2 : if (nodeOnI == List_i[ix])
381 : : {
382 : 0 : linkingEdgeNodeI = (int)ix;
383 : : }
384 [ + - ][ + - ]: 2 : else if (nodeOnI == List_ii[ix])
385 : : {
386 : 2 : linkingEdgeNodeII = (int)ix;
387 : : }
388 [ + - ][ + - ]: 2 : if (nodeOnII == List_i[ix])
389 : : {
390 : 2 : linkingEdgeNodeI = (int)ix;
391 : : }
392 [ # # ][ # # ]: 0 : else if (nodeOnII == List_ii[ix])
393 : : {
394 : 0 : linkingEdgeNodeII = (int)ix;
395 : : }
396 [ + - ][ + - ]: 2 : if (linkingEdgeNodeI != -1 && linkingEdgeNodeII != -1)
397 : : {
398 : 2 : offset = linkingEdgeNodeII - linkingEdgeNodeI;
399 [ - + ]: 2 : if (offset < 0)
400 : : {
401 : 0 : offset += size_i;
402 : : }
403 : 2 : break;
404 : : }
405 : : }
406 [ + - ][ - + ]: 2 : if (linkingEdgeNodeI == -1 || linkingEdgeNodeII == -1)
407 : : {
408 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Could not find vertices of linking surface edge on source and target.");
409 : : }
410 [ + - ][ + - ]: 2 : if (nodeOnI != List_i[linkingEdgeNodeI])
411 : : {
412 [ + - ]: 2 : std::reverse(linkEdgeNodeList.begin(), linkEdgeNodeList.end());
413 [ + - ]: 2 : nodeOnI = linkEdgeNodeList[0];
414 [ + - ]: 2 : nodeOnII = linkEdgeNodeList[linkEdgeNodeList.size() - 1];
415 : 2 : }
416 : : }
417 : : // done with all the initalizations
418 : :
419 : : // get all the position vectors in 3D
420 [ + - ][ + - ]: 4 : std::vector<Vector3D> pos_i(size_i), pos_ii(size_ii);
421 : : iGeom::Error g_err =
422 [ + - ]: 2 : mk_core()->imesh_instance()->getVtxArrCoords(&(List_i[0]),
423 [ + - ][ + - ]: 4 : size_i, iBase_INTERLEAVED, &(pos_i[0][0]));
[ + - ][ + - ]
[ + - ]
424 [ + - ]: 2 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes i");
425 [ + - ]: 2 : g_err = mk_core()->imesh_instance()->getVtxArrCoords(&(List_ii[0]),
426 [ + - ][ + - ]: 4 : size_ii, iBase_INTERLEAVED, &(pos_ii[0][0]));
[ + - ][ + - ]
[ + - ]
427 [ + - ]: 2 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes i");
428 : :
429 : : // compute the number of layers to use based on the sizing function
430 : 2 : unsigned int mesh_count = 0;
431 : : SizingFunction *surfSizing =
432 [ + - ][ + - ]: 2 : ent->mk_core()->sizing_function(ent->sizing_function_index());
[ + - ]
433 [ + - ][ - + ]: 2 : if (surfSizing->intervals() > 0)
434 : : {
435 [ # # ]: 0 : mesh_count = surfSizing->intervals();
436 [ # # ][ # # ]: 0 : if (ent->constrain_even() && mesh_count % 2)
[ # # ][ # # ]
437 : 0 : ++mesh_count;
438 : : }
439 [ + - ][ + - ]: 2 : else if (surfSizing->size() > 0)
440 : : {
441 : : double estLinkSurfWidth =
442 [ + - ][ + - ]: 2 : (boundEdges[0]->measure() + boundEdges[1]->measure()) / 2.0;
[ + - ][ + - ]
443 : : // the 1 + epsilon factor gets the correct number of edges in a test case
444 : : // where the math would work if exact, but fails due to rounding error
445 [ + - ]: 2 : double estLinkSurfLength = (1 + 1e-15) * ent->measure() / estLinkSurfWidth;
446 [ + - ]: 2 : mesh_count = estLinkSurfLength / surfSizing->size();
447 [ - + ]: 2 : if (!mesh_count)
448 : 0 : ++mesh_count;
449 [ + - ][ - + ]: 2 : if (ent->constrain_even() && mesh_count % 2)
[ # # ][ - + ]
450 : 2 : ++mesh_count;
451 : : }
452 [ # # ]: 0 : else if (linkEdgeNodeList.empty())
453 : : {
454 [ # # ]: 0 : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, "Sizing function for cylindrical linking surface with no link edge had neither positive size nor positive intervals.");
455 : : }
456 : :
457 : : // compute the interior nodes based on transforming the source and target
458 : : // edges in position
459 : 2 : unsigned int numCreatedNodes = 0;
460 [ + - ]: 2 : if (!linkEdgeNodeList.empty())
461 : : {
462 [ - + ]: 2 : if (mesh_count != linkEdgeNodeList.size() - 1)
463 : : {
464 : 0 : mesh_count = linkEdgeNodeList.size() - 1;
465 [ # # ]: 0 : std::cout << "Warning: The number of nodes on the linking surface edge "
466 [ # # ]: 0 : << "does not match the number of intervals from the sizing "
467 [ # # ]: 0 : << "function on the surface.\n";
468 : : }
469 : 2 : numCreatedNodes = (size_i - 1) * (mesh_count - 1);
470 : : }
471 : : else
472 : : {
473 : 0 : numCreatedNodes = size_i * (mesh_count - 1);
474 : : }
475 : :
476 [ + - ]: 4 : std::vector<iBase_EntityHandle> createdNodes(numCreatedNodes);
477 [ + - ]: 4 : std::vector<iBase_EntityHandle> interiorNodes(size_i * (mesh_count-1));
478 : :
479 [ + - ][ + - ]: 2 : Vector3D c0, c1;
480 [ + + ]: 4 : for (unsigned int k = 1; k < mesh_count; k++) //compute layer by layer
481 : : {
482 : 2 : double interpolationFactor = 1.0-double(k)/double(mesh_count);
483 [ + + ]: 34 : for (int i = 0; i < size_i; i++){
484 [ + - ]: 32 : c0 = pos_i[i];
485 [ + - ]: 32 : c1 = pos_ii[(i + offset) % size_i];
486 [ + - ][ + - ]: 32 : Vector3D pts = c0*interpolationFactor + c1*(1.0-interpolationFactor);
[ + - ][ + - ]
487 [ + - ]: 32 : Vector3D coords;
488 [ + - ][ + - ]: 32 : g_err = ent->igeom_instance()->getEntClosestPtTrimmed(ent->geom_handle(), pts[0], pts[1], pts[2], coords[0], coords[1], coords[2]);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
489 [ - + ]: 32 : if (g_err)
490 : : {
491 [ # # ][ # # ]: 0 : g_err = ent->igeom_instance()->getEntClosestPt(ent->geom_handle(), pts[0], pts[1], pts[2], coords[0], coords[1], coords[2]);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
492 : : }
493 [ + - ]: 32 : IBERRCHK(g_err, "Trouble get the closest xyz coordinates on the linking surface.");
494 : 32 : int pastLinkEdgeOffset = 0;
495 [ + + ]: 32 : if (i == linkingEdgeNodeI)
496 : : {
497 : : // this node corresponds to one that should already exist on the edge
498 [ + - ][ + - ]: 2 : interiorNodes[(k - 1)*size_i +i] = linkEdgeNodeList[k];
499 : : }
500 [ + - ][ + - ]: 30 : else if (i > linkingEdgeNodeI && linkingEdgeNodeI >= 0)
501 : : {
502 : 30 : pastLinkEdgeOffset = -1;
503 : : }
504 : : iMesh::Error m_err =
505 [ + - ][ + - ]: 32 : mk_core()->imesh_instance()->createVtx(coords[0], coords[1],
506 [ + - ][ + - ]: 64 : coords[2], interiorNodes[(k-1)*size_i+i]);
[ + - ][ + - ]
[ + - ]
507 [ + + ]: 32 : if (i != linkingEdgeNodeI)
508 : : {
509 [ + - ]: 30 : createdNodes[(k - 1)*(size_i - 1) + i + pastLinkEdgeOffset] =
510 [ + - ]: 30 : interiorNodes[(k - 1)*size_i + i];
511 : : }
512 [ + - ]: 32 : IBERRCHK(m_err, "Trouble create the interior node.");
513 : : }
514 : : }
515 : :
516 : : //finish creating the interior nodes
517 : : iBase_EntitySetHandle entityset;
518 [ + - ][ + - ]: 2 : iRel::Error r_err = mk_core()->irel_pair(irelPairIndex)->getEntSetRelation(ent->geom_handle(), 0, entityset);
[ + - ][ + - ]
519 [ - + ]: 2 : if (r_err){
520 : : //create the entityset
521 [ # # ][ # # ]: 0 : iMesh::Error m_err = mk_core()->imesh_instance()->createEntSet(true, entityset);
[ # # ]
522 [ # # ]: 0 : IBERRCHK(m_err, "Trouble create the entity set.");
523 : :
524 [ # # ][ # # ]: 0 : r_err = mk_core()->irel_pair(irelPairIndex)->setEntSetRelation(ent->geom_handle(), entityset);
[ # # ][ # # ]
525 [ # # ]: 0 : IBERRCHK(r_err, "Trouble create the association between the geometry and mesh entity set.");
526 : : }
527 : :
528 [ + - ][ + - ]: 2 : iMesh::Error m_err = mk_core()->imesh_instance()->addEntArrToSet(&createdNodes[0], createdNodes.size(), entityset);
[ + - ][ + - ]
529 [ + - ]: 2 : IBERRCHK(m_err, "Trouble add an array of entities to the mesh entity set.");
530 : :
531 : : // copy nodes in a vector to create quads easier
532 [ + - ]: 4 : std::vector<iBase_EntityHandle> Nodes((mesh_count+1)*size_i);
533 : : //create the int data for mesh nodes on the linking surface
534 : : iBase_TagHandle mesh_tag;
535 : : // TODO: Don't use a tag here, since multiple TFIMapping may occur at the
536 : : // same time
537 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->getTagHandle("MeshTFIMapping", mesh_tag);
[ + - ]
538 [ + + ]: 2 : if (m_err)
539 : : {
540 [ + - ][ + - ]: 1 : m_err = mk_core()->imesh_instance()->createTag("MeshTFIMapping", 1, iBase_INTEGER, mesh_tag);
[ + - ]
541 [ + - ]: 1 : IBERRCHK(m_err, "Trouble create the mesh_tag for the surface.");
542 : : }
543 : 2 : int intdata = -1;
544 [ + + ]: 34 : for (int i = 0; i < size_i; i++){
545 : 32 : intdata++;
546 [ + - ][ + - ]: 32 : m_err = mk_core()->imesh_instance()->setIntData(List_i[i], mesh_tag, intdata);
[ + - ][ + - ]
547 [ + - ]: 32 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
548 : :
549 [ + - ][ + - ]: 32 : Nodes[i] = List_i[i];
550 : :
551 [ + - ][ + - ]: 32 : m_err = mk_core()->imesh_instance()->setIntData(List_ii[i], mesh_tag, (intdata + mesh_count*size_i));
[ + - ][ + - ]
552 [ + - ]: 32 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
553 [ + - ][ + - ]: 32 : Nodes[size_i*mesh_count+i] = List_ii[i];
554 : : }
555 : :
556 [ + + ]: 34 : for (int ii = 0; ii < int(interiorNodes.size()); ii++){
557 : 32 : intdata++;
558 [ + - ][ + - ]: 32 : m_err = mk_core()->imesh_instance()->setIntData(interiorNodes[ii], mesh_tag, intdata);
[ + - ][ + - ]
559 [ + - ]: 32 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
560 : :
561 [ + - ][ + - ]: 32 : Nodes[intdata] = interiorNodes[ii];
562 : : }
563 : :
564 : : //we will always create them in the positive orientation, because we already reversed the Lists with nodes
565 [ + - ]: 4 : std::vector<iBase_EntityHandle> qNodes(4);//a generic quad
566 [ + - ]: 4 : std::vector<iBase_EntityHandle> Quads(size_i*mesh_count);
567 [ + + ]: 6 : for (unsigned int k = 0; k < mesh_count; k++){
568 [ + + ]: 68 : for (int i = 0; i < size_i; i++){
569 [ + - ][ + - ]: 64 : qNodes[0] = Nodes[ k*size_i + i ];
570 [ + - ][ + - ]: 64 : qNodes[1] = Nodes[ k*size_i + (i + 1)%size_i ];
571 [ + - ][ + - ]: 64 : qNodes[2] = Nodes[ (k+1)*size_i + (i + 1)%size_i ];
572 [ + - ][ + - ]: 64 : qNodes[3] = Nodes[ (k+1)*size_i + i ];
573 [ + - ][ + - ]: 64 : m_err = mk_core()->imesh_instance()->createEnt(iMesh_QUADRILATERAL, &qNodes[0], 4, Quads[k*size_i+i]);
[ + - ][ + - ]
[ + - ]
574 [ + - ]: 64 : IBERRCHK(m_err, "Trouble create the quadrilateral element.");
575 : : }
576 : : }
577 : :
578 : : //finish creating the quads
579 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->addEntArrToSet(&Quads[0], Quads.size(), entityset);
[ + - ][ + - ]
580 [ + - ]: 2 : IBERRCHK(m_err, "Trouble add an array of quads to the mesh entity set.");
581 : : //set int data for quads
582 [ + + ]: 66 : for (unsigned int i = 0; i < Quads.size(); i++)
583 : : {
584 [ + - ][ + - ]: 64 : m_err = mk_core()->imesh_instance()->setIntData(Quads[i], mesh_tag, i);
[ + - ][ + - ]
585 [ + - ]: 64 : IBERRCHK(m_err, "Trouble set the int data for quadrilateral elements.");
586 : : }
587 : :
588 : : //Get the global id tag
589 : 2 : const char *tag = "GLOBAL_ID";
590 : : iBase_TagHandle mesh_id_tag;
591 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->getTagHandle(tag, mesh_id_tag);
[ + - ]
592 [ + - ]: 2 : IBERRCHK(m_err, "Trouble get the mesh_id_tag for 'GLOBAL_ID'.");
593 : :
594 [ + - ][ + - ]: 4 : std::vector<iBase_EntityHandle> m_Nodes, m_Edges, m_Quads;
[ + - ]
595 : :
596 : : //set the int data for Global ID tag
597 : : iBase_EntitySetHandle root_set;
598 : : int err;
599 [ + - ][ + - ]: 2 : iMesh_getRootSet(mk_core()->imesh_instance()->instance(), &root_set, &err);
[ + - ][ + - ]
600 [ - + ]: 2 : assert(!err);
601 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_VERTEX, iMesh_POINT, m_Nodes);
[ + - ]
602 [ + - ]: 2 : IBERRCHK(m_err, "Trouble get the node list from the mesh entity set.");
603 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_EDGE, iMesh_LINE_SEGMENT, m_Edges);
[ + - ]
604 [ + - ]: 2 : IBERRCHK(m_err, "Trouble get the edges from the mesh entity set.");
605 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_FACE, iMesh_QUADRILATERAL, m_Quads);
[ + - ]
606 [ + - ]: 2 : IBERRCHK(m_err, "Trouble get the faces from the mesh entity set.");
607 : :
608 [ + + ]: 1638 : for (unsigned int i = 0; i < m_Nodes.size(); i++)
609 : : {
610 [ + - ][ + - ]: 1636 : m_err = mk_core()->imesh_instance()->setIntData(m_Nodes[i], mesh_id_tag, i);
[ + - ][ + - ]
611 [ + - ]: 1636 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
612 : : }
613 [ + + ]: 488 : for (unsigned int i = 0; i < m_Edges.size(); i++)
614 : : {
615 [ + - ][ + - ]: 486 : m_err = mk_core()->imesh_instance()->setIntData(m_Edges[i], mesh_id_tag, i);
[ + - ][ + - ]
616 [ + - ]: 486 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
617 : : }
618 [ + + ]: 1495 : for (unsigned int i = 0; i < m_Quads.size(); i++)
619 : : {
620 [ + - ][ + - ]: 1493 : m_err = mk_core()->imesh_instance()->setIntData(m_Quads[i], mesh_id_tag, i);
[ + - ][ + - ]
621 [ + - ]: 1493 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
622 : : }
623 : :
624 : : //SurfImprove(ent->geom_handle(), entityset, iBase_FACE);
625 [ + - ][ + - ]: 2 : mk_core()->save_mesh("InitialMapping.vtk");
626 : :
627 [ - + ]: 2 : if (_shapeImprove)
628 : : {
629 : : #ifdef HAVE_MESQUITE
630 : :
631 [ # # ][ # # ]: 0 : iGeom * ig_inst = mk_core()->igeom_instance(ent->iGeomIndex());
[ # # ]
632 : :
633 [ # # ][ # # ]: 0 : MeshImprove meshopt(mk_core(), true, false, false, false, ig_inst);
634 [ # # ][ # # ]: 0 : meshopt.SurfMeshImprove(ent->geom_handle(), entityset, iBase_FACE);
635 : :
636 : : #endif
637 : :
638 [ # # ][ # # ]: 0 : mk_core()->save_mesh("AfterWinslow.vtk");
639 : : #ifdef HAVE_MESQUITE
640 [ # # ][ # # ]: 0 : MeshImprove shapesmooth(mk_core(), false, false, true, false, ig_inst);
641 [ # # ][ # # ]: 0 : shapesmooth.SurfMeshImprove(ent->geom_handle(), entityset, iBase_FACE);
642 : : #endif
643 : : }
644 : :
645 : : //remove the mesh tag
646 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->rmvArrTag(&Quads[0], Quads.size(), mesh_tag);
[ + - ][ + - ]
647 [ + - ]: 2 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
648 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_i[0], List_i.size(), mesh_tag);
[ + - ][ + - ]
649 [ + - ]: 2 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
650 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_ii[0], List_ii.size(), mesh_tag);
[ + - ][ + - ]
651 [ + - ]: 2 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
652 : :
653 [ + - ]: 2 : if (interiorNodes.size() > 0)
654 : : {
655 [ + - ][ + - ]: 2 : m_err = mk_core()->imesh_instance()->rmvArrTag(&interiorNodes[0], interiorNodes.size(), mesh_tag);
[ + - ][ + - ]
656 [ + - ]: 2 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
657 : : }
658 : :
659 : 2 : return 1;
660 : : }
661 : :
662 : : /***********************************************************************************/
663 : : /*function : SurfMapping */
664 : : /*Date : Mar 3, 2011 */
665 : : /*Description: Generate the mesh on the linking surface by using TFI */
666 : : /* prepare to generate the surface by using TFI mapping interpolation */
667 : : /* 1. Get the mesh(edge mesh) from the bounding geometric edges */
668 : : /* 2. Find the corresponding relationship between edges, vertices */
669 : : /* 3. Check the nodes' corresponding relationship on the 4 bounding edges */
670 : : /* 4. Do the TFI interpolation for interior nodes' location */
671 : : /***********************************************************************************/
672 : 17 : int TFIMapping::SurfMapping(ModelEnt *ent)
673 : : {
674 : :
675 [ + - ]: 17 : int irelPairIndex = ent->iRelPairIndex();
676 [ + - ]: 17 : MEntVector boundEdges;
677 [ + - ][ + - ]: 34 : std::vector<int> senses, group_sizes;
678 [ + - ]: 17 : ent->ModelEnt::boundary(1, boundEdges, &senses, &group_sizes);
679 : : //remove this constraint in case of the side-cylinder case
680 : : //if (boundEdges.size() != 4)
681 : : // ECERRCHK(MK_FAILURE, "bad input for TFI Mapping, we can only mesh surfaces with 4 edges");
682 [ + - ]: 17 : std::cout << "Surf mapping\n";
683 : :
684 [ + - ][ + - ]: 34 : std::vector<iBase_EntityHandle> List_i, List_j, List_ii, List_jj;
[ + - ][ + - ]
685 : :
686 [ + - ]: 34 : std::vector<moab::EntityHandle> nList;
687 : : /*
688 : : corner[2] NodeList_ii -> corner[3]
689 : : ^ ^
690 : : | |
691 : : NodeList_j NodeList_jj
692 : : ^ ^
693 : : | |
694 : : corner[0] NodeList_i -> corner[1]
695 : : */
696 : : //get the nodes on each edge
697 : : // we want the start of node list i and j to be the same (corner 0)
698 [ + - ][ + - ]: 17 : boundEdges[0]->get_mesh(0, nList, true); // include start and end vertices (corners)
699 : 17 : unsigned int ix = 0;
700 : 17 : int size_i=(int)nList.size();
701 [ + + ]: 242 : for (ix = 0; ix < nList.size(); ix++)
702 [ + - ][ + - ]: 225 : List_i.push_back((iBase_EntityHandle) nList[ix]);
703 : : // if sense is reverse for edge 0, reverse list,
704 [ + - ][ + + ]: 17 : if (senses[0] == -1)
705 [ + - ]: 4 : std::reverse(List_i.begin(), List_i.end());
706 : : // so we know for sure corner 0 is at NodeList_i[0]!!
707 : 17 : nList.clear();
708 [ + - ][ + - ]: 17 : boundEdges[1]->get_mesh(0, nList, true);
709 : 17 : int size_j=(int)nList.size();
710 [ + + ]: 200 : for (ix = 0; ix < nList.size(); ix++)
711 [ + - ][ + - ]: 183 : List_jj.push_back((iBase_EntityHandle) nList[ix]);
712 [ + - ][ + + ]: 17 : if (senses[1] == -1)
713 [ + - ]: 6 : std::reverse(List_jj.begin(), List_jj.end());
714 : :
715 : 17 : nList.clear();
716 [ + - ][ + - ]: 17 : boundEdges[2]->get_mesh(0, nList, true);
717 [ + + ]: 242 : for (ix = 0; ix < nList.size(); ix++)
718 [ + - ][ + - ]: 225 : List_ii.push_back((iBase_EntityHandle) nList[ix]);
719 [ + - ][ + + ]: 17 : if (senses[2] == 1) // we reverse it if this edge is "positive" in the loop
720 [ + - ]: 9 : std::reverse(List_ii.begin(), List_ii.end());
721 : :
722 : 17 : nList.clear();
723 [ + - ][ + - ]: 17 : boundEdges[3]->get_mesh(0, nList, true);
724 [ + + ]: 200 : for (ix = 0; ix < nList.size(); ix++)
725 [ + - ][ + - ]: 183 : List_j.push_back((iBase_EntityHandle) nList[ix]);
726 [ + - ][ + + ]: 17 : if (senses[3] == 1) // we reverse it if this edge is "positive" in the loop
727 [ + - ]: 7 : std::reverse(List_j.begin(), List_j.end());
728 : :
729 [ - + ]: 17 : if (List_i.size() != List_ii.size())
730 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "opposite edges have different mesh count, abort");
731 [ - + ]: 17 : if (List_j.size() != List_jj.size())
732 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "opposite edges have different mesh count, abort");
733 : : //ok, done with all the initializations
734 : :
735 : : // get all the vectors in 3d
736 [ + - ]: 34 : std::vector<Vector3D> pos_ii(List_ii.size());
737 [ + - ]: 34 : std::vector<Vector3D> pos_i(List_i.size());
738 [ + - ]: 34 : std::vector<Vector3D> pos_j(List_j.size());
739 [ + - ]: 34 : std::vector<Vector3D> pos_jj(List_jj.size());
740 : : // iBase_INTERLEAVED
741 : : /*getVtxArrCoords( const EntityHandle* vertex_handles,
742 : : int vertex_handles_size,
743 : : StorageOrder storage_order,
744 : : double* coords_out ) const*/
745 [ + - ][ + - ]: 17 : iGeom::Error g_err = mk_core()->imesh_instance()->getVtxArrCoords(&(List_ii[0]), List_ii.size(), iBase_INTERLEAVED, &(pos_ii[0][0]));
[ + - ][ + - ]
[ + - ][ + - ]
746 [ + - ]: 17 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes ii.");
747 [ + - ][ + - ]: 17 : g_err = mk_core()->imesh_instance()->getVtxArrCoords(&(List_i[0]), List_i.size(), iBase_INTERLEAVED, &(pos_i[0][0]));
[ + - ][ + - ]
[ + - ][ + - ]
748 [ + - ]: 17 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes i.");
749 [ + - ][ + - ]: 17 : g_err = mk_core()->imesh_instance()->getVtxArrCoords(&(List_j[0]), List_j.size(), iBase_INTERLEAVED, &(pos_j[0][0]));
[ + - ][ + - ]
[ + - ][ + - ]
750 [ + - ]: 17 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes j.");
751 [ + - ][ + - ]: 17 : g_err = mk_core()->imesh_instance()->getVtxArrCoords(&(List_jj[0]), List_jj.size(), iBase_INTERLEAVED, &(pos_jj[0][0]));
[ + - ][ + - ]
[ + - ][ + - ]
752 [ + - ]: 17 : IBERRCHK(g_err, "Trouble get the xyz coordinates for nodes jj.");
753 : : //calculate the interior nodes based on transforming the top and bottom edges in position
754 [ + - ]: 34 : std::vector<iBase_EntityHandle> interiorNodes((List_j.size() - 2) * (List_i.size() - 2));
755 : : // reminder
756 : : /*
757 : : corner[2] NodeList_ii -> corner[3]
758 : : ^ ^
759 : : | |
760 : : NodeList_j NodeList_jj
761 : : ^ ^
762 : : | |
763 : : corner[0] NodeList_i -> corner[1]
764 : : */
765 [ + - ][ + - ]: 17 : Vector3D c0=pos_i[0], c1=pos_i[size_i-1], c2 = pos_ii[0], c3=pos_ii[size_i-1];
[ + - ][ + - ]
766 : :
767 [ + - ][ + - ]: 17 : Vector3D bc = 0.5*c0+0.5*c1;
[ + - ][ + - ]
768 [ + - ][ + - ]: 17 : Vector3D tc = 0.5*c2+0.5*c3;
[ + - ][ + - ]
769 [ + + ]: 166 : for (int j = 1; j < size_j - 1; j++) // compute layer by layer
770 : : // we will start from source (layer 0) to layer j (j>1, j< J-1)
771 : : // also , we will look at the target, layer J-1 to layer j
772 : : {
773 : :
774 : : // transformation from c0 and c1 to layer j
775 [ + - ][ + - ]: 149 : Matrix3D tr1 , tr2;
776 : : //target= A * ( source - 2*sc + tc) + sc
777 [ + - ][ + - ]: 149 : Vector3D cj = 0.5*pos_j[j]+0.5*pos_jj[j]; // center of layer j
[ + - ][ + - ]
[ + - ][ + - ]
778 [ + - ][ + - ]: 149 : computeTransformation(c0, c1, pos_j[j], pos_jj[j], tr1);
[ + - ]
779 : : // transformation from top, c2 and c3 to layer j
780 [ + - ][ + - ]: 149 : computeTransformation(c2, c3, pos_j[j], pos_jj[j], tr2);
[ + - ]
781 : :
782 : 149 : double interpolationFactor = j/(size_j-1.);
783 : :
784 [ + + ]: 850 : for (int i = 1; i < (size_i - 1); i++)
785 : : {
786 : : // transformation from bottom to layer j; source is b, target is j
787 [ + - ][ + - ]: 701 : Vector3D res1= tr1*(pos_i[i] -2*bc+cj)+bc;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
788 : : // transformation from top to layer j; source is t, target is j
789 [ + - ][ + - ]: 701 : Vector3D res2= tr2*(pos_ii[i] -2*tc+cj)+tc;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
790 : : // interpolate this result
791 [ + - ][ + - ]: 701 : Vector3D pts = res1*(1-interpolationFactor) + res2*interpolationFactor;
[ + - ][ + - ]
792 [ + - ]: 701 : Vector3D coords;
793 [ + - ][ + - ]: 1402 : g_err = ent->igeom_instance()->getEntClosestPtTrimmed(ent->geom_handle(), pts[0], pts[1], pts[2], coords[0], coords[1],
[ + - ][ + - ]
[ + - ][ + - ]
794 [ + - ][ + - ]: 1402 : coords[2]);
[ + - ]
795 [ + + ]: 701 : if (g_err)
796 : : {
797 [ + - ][ + - ]: 264 : g_err = ent->igeom_instance()->getEntClosestPt(ent->geom_handle(), pts[0], pts[1], pts[2], coords[0], coords[1], coords[2]);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
798 : : }
799 [ + - ]: 701 : IBERRCHK(g_err, "Trouble get the closest xyz coordinates on the linking surface.");
800 : :
801 [ + - ][ + - ]: 701 : iMesh::Error m_err = mk_core()->imesh_instance()->createVtx(coords[0], coords[1], coords[2], interiorNodes[(j - 1)
[ + - ]
802 [ + - ][ + - ]: 1402 : * (size_i - 2) + i - 1]);
[ + - ][ + - ]
803 [ + - ]: 701 : IBERRCHK(m_err, "Trouble create the interior node.");
804 : : }
805 : : }
806 : :
807 : : //finish creating the interior nodes
808 : : iBase_EntitySetHandle entityset;
809 [ + - ][ + - ]: 17 : iRel::Error r_err = mk_core()->irel_pair(irelPairIndex)->getEntSetRelation(ent->geom_handle(), 0, entityset);
[ + - ][ + - ]
810 [ - + ]: 17 : if (r_err)
811 : : {
812 : : //create the entityset
813 [ # # ][ # # ]: 0 : iMesh::Error m_err = mk_core()->imesh_instance()->createEntSet(true, entityset);
[ # # ]
814 [ # # ]: 0 : IBERRCHK(m_err, "Trouble create the entity set.");
815 : :
816 [ # # ][ # # ]: 0 : r_err = mk_core()->irel_pair(irelPairIndex)->setEntSetRelation(ent->geom_handle(), entityset);
[ # # ][ # # ]
817 [ # # ]: 0 : IBERRCHK(r_err, "Trouble create the association between the geometry and mesh entity set.");
818 : : }
819 : :
820 [ + - ][ + - ]: 17 : iMesh::Error m_err = mk_core()->imesh_instance()->addEntArrToSet(&interiorNodes[0], interiorNodes.size(), entityset);
[ + - ][ + - ]
821 [ + - ]: 17 : IBERRCHK(m_err, "Trouble add an array of entities to the mesh entity set.");
822 : :
823 : : // copy nodes in a vector to create the quads easier
824 : : // they will be arranged in layers, from bottom (j=0) towards top (j=size_j-1)
825 [ + - ]: 34 : std::vector<iBase_EntityHandle> Nodes(size_j * size_i);
826 : :
827 : : //create the int data for mesh nodes on the linking surface
828 : : iBase_TagHandle mesh_tag;
829 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->getTagHandle("MeshTFIMapping", mesh_tag);
[ + - ]
830 [ + + ]: 17 : if (m_err)
831 : : {
832 [ + - ][ + - ]: 3 : m_err = mk_core()->imesh_instance()->createTag("MeshTFIMapping", 1, iBase_INTEGER, mesh_tag);
[ + - ]
833 [ + - ]: 3 : IBERRCHK(m_err, "Trouble create the mesh_tag for the surface.");
834 : : }
835 : 17 : int intdata = -1;
836 [ + + ]: 242 : for (int i = 0; i < size_i; i++)
837 : : {
838 : 225 : intdata++;
839 [ + - ][ + - ]: 225 : m_err = mk_core()->imesh_instance()->setIntData(List_i[i], mesh_tag, intdata);
[ + - ][ + - ]
840 [ + - ]: 225 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
841 : : // bottom row, j=0
842 [ + - ][ + - ]: 225 : Nodes[i]=List_i[i];
843 : :
844 [ + - ][ + - ]: 225 : m_err = mk_core()->imesh_instance()->setIntData(List_ii[i], mesh_tag, intdata + size_i);
[ + - ][ + - ]
845 [ + - ]: 225 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
846 : : // top row, j = size_j-1
847 [ + - ][ + - ]: 225 : Nodes[ (size_i)*(size_j-1)+i] = List_ii[i];
848 : : }
849 : 17 : intdata = 2 * size_i - 1;
850 [ + + ]: 166 : for (int j = 1; j < size_j - 1; j++)
851 : : {
852 : 149 : intdata++;
853 [ + - ][ + - ]: 149 : m_err = mk_core()->imesh_instance()->setIntData(List_j[j], mesh_tag, intdata);
[ + - ][ + - ]
854 [ + - ]: 149 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
855 : : // right column, i=0
856 [ + - ][ + - ]: 149 : Nodes[size_i*j] = List_j[j];
857 : :
858 [ + - ][ + - ]: 149 : m_err = mk_core()->imesh_instance()->setIntData(List_jj[j], mesh_tag, intdata + size_j - 2);
[ + - ][ + - ]
859 [ + - ]: 149 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
860 : : // left column, i = size_i -1
861 [ + - ][ + - ]: 149 : Nodes[size_i*j + size_i-1] = List_jj[j];
862 : : }
863 : :
864 : 17 : intdata = 2 * size_i + 2 * (size_j - 2) - 1;
865 : : // it is clear that (size_i-2 > 0 and size_j-2 > 0) iff (interiorNodes.size()>0)
866 [ + + ]: 718 : for (unsigned int ii = 0; ii < interiorNodes.size(); ii++)
867 : : {
868 : 701 : intdata++;
869 [ + - ][ + - ]: 701 : m_err = mk_core()->imesh_instance()->setIntData(interiorNodes[ii], mesh_tag, intdata);
[ + - ][ + - ]
870 [ + - ]: 701 : IBERRCHK(m_err, "Trouble set the int data for mesh nodes.");
871 : 701 : int j = ii/(size_i-2) + 1;
872 : 701 : int i = (ii -(j-1)*(size_i-2)) + 1;
873 : : // copy
874 [ + - ][ + - ]: 701 : Nodes[ j*size_i + i] = interiorNodes[ii];
875 : : // compute the row and column
876 : : }
877 : :
878 : : // we will always create them in the positive orientation, because we already reversed the Lists
879 : : // with nodes
880 : :
881 [ + - ]: 34 : std::vector<iBase_EntityHandle> qNodes(4);// a generic quad
882 [ + - ]: 34 : std::vector<iBase_EntityHandle> Quads((size_j - 1) * (size_i - 1));
883 [ + + ]: 183 : for (int j=0; j <size_j-1; j++)
884 : : {
885 [ + + ]: 1224 : for (int i=0; i< size_i-1; i++)
886 : : {
887 [ + - ][ + - ]: 1058 : qNodes[0] = Nodes[ j *size_i+i ];
888 [ + - ][ + - ]: 1058 : qNodes[1] = Nodes[ j *size_i+i+1];
889 [ + - ][ + - ]: 1058 : qNodes[2] = Nodes[ (j+1)*size_i+i+1];
890 [ + - ][ + - ]: 1058 : qNodes[3] = Nodes[ (j+1)*size_i+i ];
891 [ + - ][ + - ]: 1058 : m_err = mk_core()->imesh_instance()->createEnt(iMesh_QUADRILATERAL, &qNodes[0], 4, Quads[j*(size_i-1)+i]);
[ + - ][ + - ]
[ + - ]
892 [ + - ]: 1058 : IBERRCHK(m_err, "Trouble create the quadrilateral element.");
893 : : }
894 : : }
895 : :
896 : : //finish creating the quads
897 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->addEntArrToSet(&Quads[0], Quads.size(), entityset);
[ + - ][ + - ]
898 [ + - ]: 17 : IBERRCHK(m_err, "Trouble add an array of quads to the mesh entity set.");
899 : : //set int data for quads
900 [ + + ]: 1075 : for (unsigned int i = 0; i < Quads.size(); i++)
901 : : {
902 [ + - ][ + - ]: 1058 : m_err = mk_core()->imesh_instance()->setIntData(Quads[i], mesh_tag, i);
[ + - ][ + - ]
903 [ + - ]: 1058 : IBERRCHK(m_err, "Trouble set the int data for quadrilateral elements.");
904 : : }
905 : :
906 : : //Get the global id tag
907 : 17 : const char *tag = "GLOBAL_ID";
908 : : iBase_TagHandle mesh_id_tag;
909 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->getTagHandle(tag, mesh_id_tag);
[ + - ]
910 [ + - ]: 17 : IBERRCHK(m_err, "Trouble get the mesh_id_tag for 'GLOBAL_ID'.");
911 : :
912 [ + - ][ + - ]: 34 : std::vector<iBase_EntityHandle> m_Nodes, m_Edges, m_Quads;
[ + - ]
913 : :
914 : : //set the int data for Global ID tag
915 : : iBase_EntitySetHandle root_set;
916 : : int err;
917 [ + - ][ + - ]: 17 : iMesh_getRootSet(mk_core()->imesh_instance()->instance(), &root_set, &err);
[ + - ][ + - ]
918 [ - + ]: 17 : assert(!err);
919 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_VERTEX, iMesh_POINT, m_Nodes);
[ + - ]
920 [ + - ]: 17 : IBERRCHK(m_err, "Trouble get the node list from the mesh entity set.");
921 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_EDGE, iMesh_LINE_SEGMENT, m_Edges);
[ + - ]
922 [ + - ]: 17 : IBERRCHK(m_err, "Trouble get the edges from the mesh entity set.");
923 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->getEntities(root_set, iBase_FACE, iMesh_QUADRILATERAL, m_Quads);
[ + - ]
924 [ + - ]: 17 : IBERRCHK(m_err, "Trouble get the faces from the mesh entity set.");
925 : :
926 [ + + ]: 14384 : for (unsigned int i = 0; i < m_Nodes.size(); i++)
927 : : {
928 [ + - ][ + - ]: 14367 : m_err = mk_core()->imesh_instance()->setIntData(m_Nodes[i], mesh_id_tag, i);
[ + - ][ + - ]
929 [ + - ]: 14367 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
930 : : }
931 [ + + ]: 4281 : for (unsigned int i = 0; i < m_Edges.size(); i++)
932 : : {
933 [ + - ][ + - ]: 4264 : m_err = mk_core()->imesh_instance()->setIntData(m_Edges[i], mesh_id_tag, i);
[ + - ][ + - ]
934 [ + - ]: 4264 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
935 : : }
936 [ + + ]: 11957 : for (unsigned int i = 0; i < m_Quads.size(); i++)
937 : : {
938 [ + - ][ + - ]: 11940 : m_err = mk_core()->imesh_instance()->setIntData(m_Quads[i], mesh_id_tag, i);
[ + - ][ + - ]
939 [ + - ]: 11940 : IBERRCHK(m_err, "Trouble set the int data for 'GLOBAL_ID'.");
940 : : }
941 : :
942 : : //SurfImprove(ent->geom_handle(), entityset, iBase_FACE);
943 [ + - ][ + - ]: 17 : mk_core()->save_mesh("InitialMapping.vtk");
944 : :
945 [ - + ]: 17 : if (_shapeImprove)
946 : : {
947 : : #ifdef HAVE_MESQUITE
948 : :
949 [ # # ][ # # ]: 0 : iGeom * ig_inst = mk_core()->igeom_instance(ent->iGeomIndex());
[ # # ]
950 : :
951 [ # # ][ # # ]: 0 : MeshImprove meshopt(mk_core(), true, false, false, false, ig_inst);
952 [ # # ][ # # ]: 0 : meshopt.SurfMeshImprove(ent->geom_handle(), entityset, iBase_FACE);
953 : :
954 : : #endif
955 : : //mk_core()->save_mesh("AfterLaplace.vtk");
956 : :
957 : : //if there is the parametric space, let Winslow smooth inside the parametric space
958 [ # # ]: 0 : SmoothWinslow(List_i, List_ii, List_j, List_jj, interiorNodes, Quads, mesh_tag, ent);
959 : :
960 [ # # ][ # # ]: 0 : mk_core()->save_mesh("AfterWinslow.vtk");
961 : : #ifdef HAVE_MESQUITE
962 [ # # ][ # # ]: 0 : MeshImprove shapesmooth(mk_core(), false, false, true, false, ig_inst);
963 [ # # ][ # # ]: 0 : shapesmooth.SurfMeshImprove(ent->geom_handle(), entityset, iBase_FACE);
964 : : #endif
965 : : }
966 : : //remove the mesh tag
967 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&Quads[0], Quads.size(), mesh_tag);
[ + - ][ + - ]
968 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
969 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_i[0], List_i.size(), mesh_tag);
[ + - ][ + - ]
970 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
971 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_ii[0], List_ii.size(), mesh_tag);
[ + - ][ + - ]
972 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
973 [ + - ]: 17 : if (List_j.size() > 2)
974 : : {
975 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_j[0], List_j.size() - 2, mesh_tag);
[ + - ][ + - ]
976 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
977 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&List_jj[0], List_jj.size() - 2, mesh_tag);
[ + - ][ + - ]
978 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
979 : : }
980 [ + - ]: 17 : if (interiorNodes.size() > 0)
981 : : {
982 [ + - ][ + - ]: 17 : m_err = mk_core()->imesh_instance()->rmvArrTag(&interiorNodes[0], interiorNodes.size(), mesh_tag);
[ + - ][ + - ]
983 [ + - ]: 17 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
984 : : }
985 : :
986 : 17 : return 1;
987 : : }
988 : 298 : void TFIMapping::computeTransformation(Vector3D & A, Vector3D & B, Vector3D & C, Vector3D & D,
989 : : Matrix3D & M)
990 : : {
991 [ + - ]: 298 : Matrix3D tmpMatrix;
992 [ + - ]: 298 : Matrix3D bMatrix;
993 [ + - ][ + - ]: 298 : Vector3D c1=0.5*A+0.5*B;
[ + - ][ + - ]
994 [ + - ][ + - ]: 298 : Vector3D c2=0.5*C+0.5*D;
[ + - ][ + - ]
995 [ + - ][ + - ]: 298 : Vector3D normal=(A-D)*(B-C);// this should be not modified by the transformation
[ + - ][ + - ]
996 : : // we add this just to increase the rank of tmpMatrix
997 : : // the normal to the "plane" of interest should not be rotated at all
998 : : // as if we say that we want A*normal = normal
999 [ + - ][ + - ]: 298 : Vector3D s1=A-2*c1+c2;
[ + - ][ + - ]
1000 [ + - ][ + - ]: 298 : Vector3D s2=B-2*c1+c2;
[ + - ][ + - ]
1001 [ + - ][ + - ]: 298 : Vector3D t1=C-c1;
1002 [ + - ][ + - ]: 298 : Vector3D t2=D-c1;
1003 : : // so we are looking for M such that
1004 : : /*
1005 : : * M*s1 = t1
1006 : : * M*s2 = t2
1007 : : * M*n = n
1008 : : *
1009 : : * In simple cases, M is identity
1010 : : */
1011 : :
1012 [ + - ]: 298 : tmpMatrix.set_column(0, s1);
1013 [ + - ]: 298 : tmpMatrix.set_column(1, s2);
1014 [ + - ]: 298 : tmpMatrix.set_column(2, normal);
1015 [ + - ]: 298 : bMatrix.set_column(0, t1);
1016 [ + - ]: 298 : bMatrix.set_column(1, t2);
1017 [ + - ]: 298 : bMatrix.set_column(2, normal);
1018 : :
1019 [ + - ]: 298 : double detValue = det(tmpMatrix);
1020 : : (void) detValue;
1021 [ - + ]: 298 : assert(detValue*detValue>1.e-20);
1022 : :
1023 : : //solve the affine mapping matrix, make use of inverse matrix to get affine mapping matrix
1024 [ + - ]: 298 : Matrix3D InvMatrix = inverse(tmpMatrix);
1025 [ + - ]: 298 : M = bMatrix*InvMatrix;
1026 : :
1027 : 298 : }
1028 : : //smooth the quadrilateral mesh on the linking surface
1029 : 0 : void TFIMapping::SmoothWinslow(std::vector<iBase_EntityHandle> &List_i, std::vector<iBase_EntityHandle> &List_ii, std::vector<
1030 : : iBase_EntityHandle> &List_j, std::vector<iBase_EntityHandle> &List_jj, std::vector<iBase_EntityHandle> &interiorNodes,
1031 : : std::vector<iBase_EntityHandle> &quads, iBase_TagHandle &taghandle, ModelEnt *ent)
1032 : : {
1033 [ # # ]: 0 : std::vector<std::set<int> > AdjElements;
1034 [ # # ]: 0 : std::vector<std::vector<int> > Quads;
1035 [ # # ]: 0 : std::vector<std::vector<double> > coords;
1036 [ # # ]: 0 : std::vector<bool> isBnd;
1037 [ # # ]: 0 : std::vector<iBase_EntityHandle> nodes;
1038 [ # # ]: 0 : std::vector<double> weight;
1039 : :
1040 : 0 : bool isParameterized = false;
1041 : :
1042 [ # # ][ # # ]: 0 : iGeom::Error g_err = ent->igeom_instance()->isEntParametric(ent->geom_handle(), isParameterized);
[ # # ]
1043 [ # # ]: 0 : IBERRCHK(g_err, "Trouble check whether the surface is parameterized or not.");
1044 : 0 : isParameterized = false;
1045 : :
1046 : : //resize the coords to store all the nodes's coordinates on the linking surface
1047 [ # # ]: 0 : coords.resize(List_i.size() * List_j.size());
1048 [ # # ]: 0 : isBnd.resize(coords.size());
1049 [ # # ]: 0 : nodes.resize(coords.size());
1050 [ # # ]: 0 : for (unsigned int i = 0; i < coords.size(); i++)
1051 [ # # ][ # # ]: 0 : coords[i].resize(3);
1052 : :
1053 : : iMesh::Error m_err;
1054 : : //input the boundary nodes
1055 [ # # ]: 0 : for (unsigned int i = 0; i < List_i.size(); i++)
1056 : : {
1057 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getVtxCoord(List_i[i], coords[i][0], coords[i][1], coords[i][2]);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1058 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the vertex coordinates.");
1059 [ # # ][ # # ]: 0 : nodes[i] = List_i[i];
1060 [ # # ]: 0 : isBnd[i] = true;
1061 : :
1062 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getVtxCoord(List_ii[i], coords[List_i.size() + i][0], coords[List_i.size() + i][1],
[ # # ][ # # ]
[ # # ]
1063 [ # # ][ # # ]: 0 : coords[List_i.size() + i][2]);
[ # # ][ # # ]
[ # # ]
1064 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the vertex coordinates.");
1065 [ # # ]: 0 : if (isParameterized)
1066 : : {
1067 : : double uv[2];
1068 [ # # ][ # # ]: 0 : g_err = ent->igeom_instance()->getEntXYZtoUV(ent->geom_handle(), coords[List_i.size() + i][0], coords[List_i.size() + i][1],
[ # # ][ # # ]
[ # # ]
1069 [ # # ][ # # ]: 0 : coords[List_i.size() + i][2], uv[0], uv[1]);
[ # # ][ # # ]
1070 [ # # ]: 0 : IBERRCHK(g_err, "Trouble get the uv from xyz.");
1071 [ # # ][ # # ]: 0 : coords[List_i.size() + i][0] = uv[0];
1072 [ # # ][ # # ]: 0 : coords[List_i.size() + i][1] = uv[1];
1073 : : }
1074 : :
1075 [ # # ][ # # ]: 0 : nodes[List_i.size() + i] = List_ii[i];
1076 [ # # ]: 0 : isBnd[List_i.size() + i] = true;
1077 : : }
1078 [ # # ]: 0 : if (int(List_j.size()) > 2)
1079 : : {
1080 [ # # ]: 0 : for (unsigned int i = 1; i < (List_j.size() - 1); i++)
1081 : : {
1082 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getVtxCoord(List_j[i], coords[2 * List_i.size() + i - 1][0], coords[2 * List_i.size()
[ # # ]
1083 [ # # ]: 0 : + i - 1][1], coords[2 * List_i.size() + i - 1][2]);
[ # # # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1084 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the vertex coordinates.");
1085 [ # # ][ # # ]: 0 : nodes[2 * List_i.size() + i - 1] = List_j[i];
1086 [ # # ]: 0 : isBnd[2 * List_i.size() + i - 1] = true;
1087 : :
1088 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getVtxCoord(List_jj[i], coords[2 * List_i.size() + List_j.size() - 2 + i - 1][0],
[ # # ]
1089 [ # # ][ # # ]: 0 : coords[2 * List_i.size() + List_j.size() - 2 + i - 1][1], coords[2 * List_i.size() + List_j.size() - 2 + i - 1][2]);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1090 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the vertex coordinates.");
1091 [ # # ][ # # ]: 0 : nodes[2 * List_i.size() + List_j.size() - 2 + i - 1] = List_jj[i];
1092 [ # # ]: 0 : isBnd[2 * List_i.size() + List_j.size() - 2 + i - 1] = true;
1093 : : }
1094 : : }
1095 : : //input the interior nodes
1096 [ # # ]: 0 : if (interiorNodes.size() > 0)
1097 : : {
1098 [ # # ]: 0 : for (unsigned int i = 0; i < interiorNodes.size(); i++)
1099 : : {
1100 [ # # ]: 0 : m_err = mk_core()->imesh_instance()->getVtxCoord(interiorNodes[i],
1101 [ # # ][ # # ]: 0 : coords[2 * List_i.size() + 2 * (List_j.size() - 2) + i][0], coords[2 * List_i.size() + 2 * (List_j.size() - 2) + i][1],
[ # # ][ # # ]
1102 [ # # ][ # # ]: 0 : coords[2 * List_i.size() + 2 * (List_j.size() - 2) + i][2]);
[ # # ][ # # ]
[ # # ]
1103 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the vertex coordinates.");
1104 [ # # ][ # # ]: 0 : nodes[2 * List_i.size() + 2 * (List_j.size() - 2) + i] = interiorNodes[i];
1105 [ # # ]: 0 : isBnd[2 * List_i.size() + 2 * (List_j.size() - 2) + i] = false;
1106 : : }
1107 : : }
1108 : :
1109 : : //update the AdjElements info
1110 : : //notice: during this process, adjacent quads will be returned around a node. The quads from source surface and target surface may be returned.
1111 [ # # ]: 0 : AdjElements.resize(nodes.size());
1112 [ # # ]: 0 : for (unsigned int i = 0; i < AdjElements.size(); i++)
1113 : : {
1114 [ # # ][ # # ]: 0 : if (!isBnd[i])
1115 : : {
1116 [ # # ]: 0 : std::vector<iBase_EntityHandle> adjEnts;
1117 : 0 : adjEnts.clear();
1118 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getEntAdj(nodes[i], iBase_FACE, adjEnts);
[ # # ][ # # ]
1119 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the adjacent quads wrt a node.");
1120 [ # # ]: 0 : for (unsigned int j = 0; j < adjEnts.size(); j++)
1121 : : {
1122 : 0 : int index_id = -1;
1123 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getIntData(adjEnts[j], taghandle, index_id);
[ # # ][ # # ]
1124 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get int data for quads.");
1125 [ # # ][ # # ]: 0 : AdjElements[i].insert(index_id);
1126 : : //std::cout<< " i= " << i << " index_id:" << index_id << "\n";
1127 : 0 : }
1128 : : }
1129 : : }
1130 : :
1131 : : //update the Quads' info
1132 [ # # ]: 0 : Quads.resize(quads.size());
1133 [ # # ]: 0 : for (unsigned int i = 0; i < Quads.size(); i++)
1134 : : {
1135 [ # # ]: 0 : std::vector<iBase_EntityHandle> adjEnts;
1136 : 0 : adjEnts.clear();
1137 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getEntAdj(quads[i], iBase_VERTEX, adjEnts);
[ # # ][ # # ]
1138 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the adjacent nodes wrt a quad.");
1139 : :
1140 [ # # ]: 0 : assert(adjEnts.size()==4);
1141 [ # # ][ # # ]: 0 : Quads[i].resize(4);
1142 : :
1143 [ # # ]: 0 : for (unsigned int j = 0; j < adjEnts.size(); j++)
1144 : : {
1145 : 0 : int index_id = -1;
1146 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getIntData(adjEnts[j], taghandle, index_id);
[ # # ][ # # ]
1147 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get int data for nodes.");
1148 [ # # ][ # # ]: 0 : Quads[i][j] = index_id;
1149 : : }
1150 : 0 : }
1151 : :
1152 : : //detect the connectivity
1153 [ # # ][ # # ]: 0 : std::vector<std::vector<int> > connect(nodes.size(), std::vector<int>(9));
1154 [ # # ]: 0 : for (unsigned int i = 0; i < nodes.size(); i++)
1155 : : {
1156 [ # # ][ # # ]: 0 : if (!isBnd[i])
1157 : : {
1158 : : //there are 4 adjacent quadrilateral elements around node i
1159 : : //std::cout << " element i:" << i << " AdjElements[i].size() " << AdjElements[i].size() << "\n";
1160 [ # # ][ # # ]: 0 : assert(AdjElements[i].size() == 4);
1161 [ # # ]: 0 : std::set<int>::iterator it = AdjElements[i].begin();
1162 : : int st_index[4];
1163 : : //process 4 quad elements
1164 : 0 : int j = -1;
1165 [ # # ][ # # ]: 0 : for (; it != AdjElements[i].end(); it++)
[ # # ][ # # ]
1166 : : {
1167 : 0 : j++;
1168 [ # # ][ # # ]: 0 : if (int(i) == Quads[*it][0])
[ # # ][ # # ]
1169 : 0 : st_index[j] = 0;
1170 [ # # ][ # # ]: 0 : else if (int(i) == Quads[*it][1])
[ # # ][ # # ]
1171 : 0 : st_index[j] = 1;
1172 [ # # ][ # # ]: 0 : else if (int(i) == Quads[*it][2])
[ # # ][ # # ]
1173 : 0 : st_index[j] = 2;
1174 : : else
1175 : 0 : st_index[j] = 3;
1176 : : }
1177 [ # # ]: 0 : it = AdjElements[i].begin();
1178 [ # # ][ # # ]: 0 : connect[i][2] = Quads[*it][(st_index[0] + 3) % 4];
[ # # ][ # # ]
[ # # ]
1179 [ # # ][ # # ]: 0 : connect[i][8] = Quads[*it][(st_index[0] + 1) % 4];
[ # # ][ # # ]
[ # # ]
1180 [ # # ][ # # ]: 0 : connect[i][1] = Quads[*it][(st_index[0] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1181 : : //finish processing the quad 1
1182 [ # # ]: 0 : std::set<int>::iterator it1 = AdjElements[i].begin();
1183 [ # # ]: 0 : it1++;
1184 [ # # ][ # # ]: 0 : for (j = 1; j < 4; j++, it1++)
1185 : : {
1186 [ # # ][ # # ]: 0 : if (connect[i][8] == Quads[*it1][(st_index[j] + 1) % 4])
[ # # ][ # # ]
[ # # ][ # # ]
1187 : : {
1188 [ # # ][ # # ]: 0 : connect[i][7] = Quads[*it1][(st_index[j] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1189 [ # # ][ # # ]: 0 : connect[i][6] = Quads[*it1][(st_index[j] + 3) % 4];
[ # # ][ # # ]
[ # # ]
1190 : 0 : break;
1191 : : }
1192 [ # # ][ # # ]: 0 : else if (connect[i][8] == Quads[*it1][(st_index[j] + 3) % 4])
[ # # ][ # # ]
[ # # ][ # # ]
1193 : : {
1194 [ # # ][ # # ]: 0 : connect[i][7] = Quads[*it1][(st_index[j] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1195 [ # # ][ # # ]: 0 : connect[i][6] = Quads[*it1][(st_index[j] + 1) % 4];
[ # # ][ # # ]
[ # # ]
1196 : 0 : break;
1197 : : }
1198 : : else
1199 : 0 : continue;
1200 : : }
1201 : : //finish processing the quad 2
1202 [ # # ]: 0 : std::set<int>::iterator it2 = AdjElements[i].begin();
1203 [ # # ]: 0 : it2++;
1204 [ # # ][ # # ]: 0 : for (j = 1; it2 != AdjElements[i].end(); it2++, j++)
[ # # ][ # # ]
1205 : : {
1206 [ # # ][ # # ]: 0 : if (connect[i][2] == Quads[*it2][(st_index[j] + 1) % 4])
[ # # ][ # # ]
[ # # ][ # # ]
1207 : : {
1208 [ # # ][ # # ]: 0 : connect[i][3] = Quads[*it2][(st_index[j] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1209 [ # # ][ # # ]: 0 : connect[i][4] = Quads[*it2][(st_index[j] + 3) % 4];
[ # # ][ # # ]
[ # # ]
1210 : 0 : break;
1211 : : }
1212 [ # # ][ # # ]: 0 : else if (connect[i][2] == Quads[*it2][(st_index[j] + 3) % 4])
[ # # ][ # # ]
[ # # ][ # # ]
1213 : : {
1214 [ # # ][ # # ]: 0 : connect[i][3] = Quads[*it2][(st_index[j] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1215 [ # # ][ # # ]: 0 : connect[i][4] = Quads[*it2][(st_index[j] + 1) % 4];
[ # # ][ # # ]
[ # # ]
1216 : 0 : break;
1217 : : }
1218 : : else
1219 : 0 : continue;
1220 : : }
1221 : : //finish processing the quad 4;
1222 [ # # ]: 0 : std::set<int>::iterator it3 = AdjElements[i].begin();
1223 [ # # ]: 0 : it3++;
1224 [ # # ][ # # ]: 0 : for (j = 1; it3 != AdjElements[i].end(); it3++, j++)
[ # # ][ # # ]
1225 : : {
1226 [ # # ][ # # ]: 0 : if ((it3 != it1) && (it3 != it2))
[ # # ][ # # ]
[ # # ]
1227 : : {
1228 [ # # ][ # # ]: 0 : connect[i][5] = Quads[*it2][(st_index[j] + 2) % 4];
[ # # ][ # # ]
[ # # ]
1229 : : }
1230 : : else
1231 : 0 : continue;
1232 : : }
1233 : : }
1234 : : }
1235 : : //finish all the initialization
1236 : :
1237 [ # # ]: 0 : EquipotentialSmooth smoother;
1238 : :
1239 : : //IsoLaplace smoother;
1240 : :
1241 [ # # ][ # # ]: 0 : smoother.SetupData(AdjElements, Quads, coords, isBnd, connect);
[ # # ][ # # ]
[ # # ][ # # ]
1242 [ # # ]: 0 : smoother.Execute();
1243 : :
1244 : : //std::vector<std::vector<double> > coors;
1245 [ # # ]: 0 : smoother.GetCoords(coords);
1246 : :
1247 : : //update the new position for nodes
1248 [ # # ]: 0 : for (unsigned int i = 0; i < nodes.size(); i++)
1249 : : {
1250 [ # # ][ # # ]: 0 : if (!isBnd[i])
1251 : : {
1252 [ # # ][ # # ]: 0 : double tmp_coord[3] = { coords[i][0], coords[i][1], coords[i][2] };
[ # # ][ # # ]
[ # # ][ # # ]
1253 [ # # ]: 0 : if (!isParameterized)
1254 : : {
1255 [ # # ][ # # ]: 0 : iGeom::Error g_err = ent->igeom_instance()->getEntClosestPt(ent->geom_handle(), coords[i][0], coords[i][1], coords[i][2],
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1256 [ # # ][ # # ]: 0 : tmp_coord[0], tmp_coord[1], tmp_coord[2]);
1257 [ # # ]: 0 : IBERRCHK(g_err, "Trouble get a closest point on the linking surface.");
1258 : : }
1259 : : else
1260 : : {
1261 [ # # ][ # # ]: 0 : iGeom::Error g_err = ent->igeom_instance()->getEntXYZtoUV(ent->geom_handle(), coords[i][0], coords[i][1], tmp_coord[0],
[ # # ][ # # ]
[ # # ]
1262 [ # # ][ # # ]: 0 : tmp_coord[1], tmp_coord[2]);
1263 [ # # ]: 0 : IBERRCHK(g_err, "Trouble get the xyz from uv.");
1264 : :
1265 : : }
1266 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->setVtxCoord(nodes[i], tmp_coord[0], tmp_coord[1], tmp_coord[2]);
[ # # ][ # # ]
1267 [ # # ]: 0 : IBERRCHK(m_err, "Trouble set the new coordinates for nodes.");
1268 : : }
1269 : : }
1270 : :
1271 : : //remove the unnecessary tag after smoothing
1272 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->rmvArrTag(&nodes[0], nodes.size(), taghandle);
[ # # ][ # # ]
1273 [ # # ]: 0 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
1274 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->rmvArrTag(&quads[0], quads.size(), taghandle);
[ # # ][ # # ]
1275 [ # # ]: 0 : IBERRCHK(m_err, "Trouble remove the tag values from an array of entities.");
1276 : : //m_err = mk_core()->imesh_instance()->destroyTag(taghandle, 1);
1277 : : //IBERRCHK(m_err, "Trouble destroy a tag.");
1278 : 0 : }
1279 : :
1280 [ + - ][ + - ]: 156 : } // namespace MeshKit
1281 : :
|
