Branch data Line data Source code
1 : : #ifdef _WIN32
2 : : #define _USE_MATH_DEFINES
3 : : #endif
4 : :
5 : : #include "meshkit/OneToOneSwept.hpp"
6 : : #include "meshkit/iMesh.hpp"
7 : : #include "meshkit/EdgeMesher.hpp"
8 : : #include "meshkit/TFIMapping.hpp"
9 : : #ifdef HAVE_MESQUITE
10 : : #include "meshkit/MeshImprove.hpp"
11 : : #endif
12 : : #include <iostream>
13 : : #include <algorithm>
14 : : #include <math.h>
15 : : #include <map>
16 : : #ifdef HAVE_ARMADILLO
17 : : #include "SurfProHarmonicMap.hpp"
18 : : #endif
19 : :
20 : : namespace MeshKit {
21 : :
22 : : //---------------------------------------------------------------------------//
23 : : //Entity Type initialization for OneToOneSwept meshing
24 : : moab::EntityType OneToOneSwept_tps[] = { moab::MBVERTEX, moab::MBQUAD, moab::MBHEX, moab::MBMAXTYPE };
25 : 40 : const moab::EntityType* OneToOneSwept::output_types()
26 : : {
27 : 40 : return OneToOneSwept_tps;
28 : : }
29 : :
30 : : //---------------------------------------------------------------------------//
31 : : // construction function for OneToOneSwept class
32 : 6 : OneToOneSwept::OneToOneSwept(MKCore *mk_core, const MEntVector &me_vec) :
33 [ + - ][ + - ]: 6 : MeshScheme(mk_core, me_vec)
[ + - ][ + - ]
34 : : {
35 : 6 : }
36 : :
37 : : //---------------------------------------------------------------------------//
38 : : // deconstruction function for OneToOneSwept class
39 : 18 : OneToOneSwept::~OneToOneSwept()
40 : : {
41 [ - + ]: 12 : }
42 : :
43 : : // these have to be called before setup, otherwise we don't know source and target
44 : : //---------------------------------------------------------------------------//
45 : : // set the source surface function
46 : 6 : void OneToOneSwept::SetSourceSurface(int index)
47 : : {
48 : 6 : index_src = index;
49 : 6 : }
50 : :
51 : : //---------------------------------------------------------------------------//
52 : : // set the target surface function
53 : 6 : void OneToOneSwept::SetTargetSurface(int index)
54 : : {
55 : 6 : index_tar = index;
56 : 6 : }
57 : :
58 : : //---------------------------------------------------------------------------//
59 : : // setup function: define the size between the different layers
60 : 6 : void OneToOneSwept::setup_this()
61 : : {
62 : : //compute the number of intervals for the associated ModelEnts, from the size set on them
63 : : //the sizing function they point to, or a default sizing function
64 : :
65 [ - + ]: 6 : if (mentSelection.empty())
66 : 6 : return;
67 [ + - ]: 6 : ModelEnt * firstME = (mentSelection.begin())->first;
68 : 6 : igeom_inst = firstME->igeom_instance();
69 : : /*moab::Interface **/
70 : 6 : mb = mk_core()->moab_instance();// we should get it also from
71 : : // model entity, if we have multiple moab instances running around
72 : :
73 : 6 : const char *tag = "GLOBAL_ID";
74 : 6 : iGeom::Error g_err = igeom_inst->getTagHandle(tag, geom_id_tag);
75 : 6 : IBERRCHK(g_err, "Trouble get the geom_id_tag for 'GLOBAL_ID'.");
76 : :
77 [ + - ][ + - ]: 12 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
78 : : {
79 [ + - ]: 6 : ModelEnt *me = mit->first;
80 : : //first check to see whether entity is meshed
81 [ + - ][ + - ]: 6 : if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0)
[ + - ][ - + ]
[ - + ]
82 : 0 : continue;
83 : :
84 [ + - ][ + - ]: 6 : SizingFunction *sf = mk_core()->sizing_function(me->sizing_function_index());
[ + - ]
85 [ - + ][ # # ]: 6 : if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT && mk_core()->sizing_function(0))
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ - + ]
86 [ # # ][ # # ]: 0 : sf = mk_core()->sizing_function(0);
87 : :
88 [ - + ][ # # ]: 6 : if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT)
[ # # ][ # # ]
[ # # ][ - + ]
89 : : {
90 : : //no sizing set, just assume default #intervals as 4
91 [ # # ]: 0 : me->mesh_intervals(4);
92 [ # # ]: 0 : me->interval_firmness(DEFAULT);
93 : : }
94 : : else
95 : : {
96 : : //check # intervals first, then size, and just choose for now
97 [ + - ][ + + ]: 6 : if (sf->intervals() > 0)
98 : : {
99 [ + - ][ - + ]: 2 : if (me->constrain_even() && sf->intervals() % 2)
[ # # ][ # # ]
[ - + ]
100 [ # # ][ # # ]: 0 : me -> mesh_intervals(sf->intervals() + 1);
101 : : else
102 [ + - ][ + - ]: 2 : me -> mesh_intervals(sf->intervals());
103 [ + - ]: 2 : me -> interval_firmness(HARD);
104 : : }
105 [ + - ][ + - ]: 4 : else if (sf->size() > 0)
106 : : {
107 : : // FIXME: This calculation is bogus if size means edge lengths.
108 : : // This should divide by the measurement of the source surface.
109 [ + - ][ + - ]: 4 : int intervals = me->measure() / sf->size();
110 [ - + ]: 4 : if (!intervals)
111 : 0 : intervals++;
112 [ + - ][ - + ]: 4 : if (me->constrain_even() && intervals % 2)
[ # # ][ - + ]
113 : 0 : intervals++;
114 [ + - ]: 4 : me->mesh_intervals(intervals);
115 [ + - ]: 4 : me->interval_firmness(SOFT);
116 : : }
117 : : else
118 [ # # ]: 0 : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, "Sizing function for edge had neither positive size nor positive intervals.");
119 : : }
120 : :
121 [ + - ]: 6 : int numIntervals = me->mesh_intervals();
122 [ + - ][ + - ]: 6 : if (!sf || sf->intervals() != numIntervals)
[ + + ][ + + ]
123 : : {
124 [ + - ][ + - ]: 4 : sf = new SizingFunction(mk_core(), numIntervals, -1.);
[ + - ]
125 : : }
126 : :
127 [ + - ]: 6 : std::vector<iBase_EntityHandle> gFaces;
128 [ + - ][ + - ]: 6 : g_err = igeom_inst->getEntAdj(me->geom_handle(), iBase_FACE, gFaces);
129 [ + - ]: 6 : IBERRCHK(g_err, "Trouble get the geometrical adjacent to the solid");
130 : : //select the source surface and target surface
131 : : // these were setup from the beginning, it could be different for each volume!!!!
132 : : // we need a better way to decide/select target and source surfaces
133 : : // this also assumes that the sourceSurface is fixed between setup and execute
134 [ + - ]: 6 : sourceSurface = gFaces[index_src];
135 [ + - ]: 6 : targetSurface = gFaces[index_tar];
136 [ + - ]: 12 : MEntVector linkingSurfaces;
137 : : int index_id_src, index_id_tar;
138 [ + - ]: 6 : iGeom::Error g_err = igeom_inst->getIntData(sourceSurface, geom_id_tag, index_id_src);
139 [ + - ]: 6 : IBERRCHK(g_err, "Trouble get the int data for source surface.");
140 [ + - ]: 6 : g_err = igeom_inst->getIntData(targetSurface, geom_id_tag, index_id_tar);
141 [ + - ]: 6 : IBERRCHK(g_err, "Trouble get the int data for target surface.");
142 [ + - ][ + - ]: 6 : std::cout << " source Global ID: " << index_id_src << " target Global ID: " << index_id_tar << "\n";
[ + - ][ + - ]
[ + - ]
143 : :
144 [ + + ]: 42 : for (unsigned int i = 0; i < gFaces.size(); i++)
145 : : {
146 : : int gid;
147 [ + - ][ + - ]: 36 : g_err = igeom_inst->getIntData(gFaces[i], geom_id_tag, gid);
148 [ + - ]: 36 : IBERRCHK(g_err, "Trouble get the int data for source surface.");
149 [ + - ]: 36 : std::vector<iBase_EntityHandle> gEdges;
150 [ + - ][ + - ]: 36 : g_err = igeom_inst->getEntAdj(gFaces[i], iBase_EDGE, gEdges);
151 [ + - ]: 36 : IBERRCHK(g_err, "Trouble get the geometrical edges adjacent to the surface");
152 [ + - ][ + - ]: 36 : std::cout << " face index " << i << " with ID: " << gid << " has " << gEdges.size() << " edges\n";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
153 : 36 : }
154 [ + - ]: 12 : MEntVector surfs;
155 [ + - ]: 6 : me->get_adjacencies(2, surfs);// all surfaces adjacent to the volume
156 [ + + ]: 42 : for (unsigned int i = 0; i < surfs.size(); i++)
157 : : {
158 : : int index_id_link;
159 [ + - ][ + - ]: 36 : g_err = igeom_inst->getIntData(surfs[i]->geom_handle(), geom_id_tag, index_id_link);
[ + - ]
160 [ + - ]: 36 : IBERRCHK(g_err, "Trouble get the int data for linking surface.");
161 [ + + ][ + + ]: 36 : if ((index_id_link != index_id_src) && (index_id_link != index_id_tar))
162 : : {
163 [ + - ]: 24 : MEntVector edges;
164 [ + - ][ + - ]: 24 : surfs[i]->get_adjacencies(1, edges);// all surfaces adjacent to the volume
165 [ + - ][ + - ]: 24 : std::cout << " linking surface with global ID: " << index_id_link << " has " << edges.size() << " edges\n";
[ + - ][ + - ]
[ + - ]
166 : : //assert((int)edges.size()==4);
167 : :
168 [ + - ][ + - ]: 24 : linkingSurfaces.push_back(surfs[i]);
169 : : }
170 : : }
171 : : // now create a TFI mapping
172 : : // this will trigger edge meshers for linking edges, and for target surface edges
173 [ + - ][ + - ]: 6 : TFIMapping *tm = (TFIMapping*) mk_core()->construct_meshop("TFIMapping", linkingSurfaces);
[ + - ]
174 : : // For now, do not assign any sizing function at all
175 : : // for (unsigned int i = 0; i < linkingSurfaces.size(); i++)
176 : : // {
177 : : // linkingSurfaces[i]->sizing_function_index(sf->core_index());
178 : : // }
179 [ + - ][ + - ]: 6 : mk_core()->insert_node(tm, (MeshOp*) this);
180 : 6 : }
181 : :
182 : 6 : mk_core()->print_graph("AfterOneSetup.eps");
183 : :
184 : : }
185 : :
186 : : //---------------------------------------------------------------------------//
187 : : // execute function: generate the all-hex mesh through sweeping from source
188 : : // surface to target surface
189 : 6 : void OneToOneSwept::execute_this()
190 : : {
191 : :
192 [ + - ][ + - ]: 12 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
193 : : {
194 [ + - ]: 6 : ModelEnt *me = mit -> first;
195 : : // maybe it will be decided later?
196 : : //case 1: if numLayers = 1, then it is not necessary to create any vertices for linking surface, All the vertices have been created by source surface and target surface
197 : : // bLayers will contain nodes in order, layer by layer
198 [ + - ]: 6 : std::vector<moab::EntityHandle> bLayers;
199 [ + - ]: 6 : BuildLateralBoundaryLayers(me, bLayers);
200 : : //int sizeBLayer = bLayers.size()/(numLayers+1);
201 [ + - ]: 6 : volume = me->geom_handle();
202 : :
203 [ + - ]: 6 : TargetSurfProjection(bLayers);// the top layer is the last list of nodes
204 : :
205 : : //get the volume mesh entity set
206 [ + - ][ + - ]: 6 : iRel::Error r_err = mk_core()->irel_pair(me->iRelPairIndex())->getEntSetRelation(me->geom_handle(), 0, volumeSet);
[ + - ][ + - ]
[ + - ]
207 : 6 : IBERRCHK(r_err,
208 [ + - ]: 6 : "Trouble get the volume mesh entity set from the geometrical volume.");
209 : :
210 [ + - ][ + - ]: 12 : vector<vector<Vertex> > linkVertexList(numLayers - 1, vector<Vertex> (NodeList.size()));
211 : : // this will compute the interior points distribution, in each layer
212 [ + - ]: 6 : ProjectInteriorLayers(bLayers, linkVertexList);
213 [ + - ][ + - ]: 6 : std::cout << "Projected interior layers." << std::endl;
214 : :
215 [ + - ]: 6 : CreateElements(linkVertexList);
216 [ + - ][ + - ]: 6 : std::cout << "Created elements." << std::endl;
217 : :
218 [ + - ][ + - ]: 6 : mk_core()->save_mesh("BeforeVolumeImprove.h5m");
219 : : #if HAVE_MESQUITE
220 : : //MeshImprove meshImpr(mk_core());
221 : : //meshImpr.VolumeMeshImprove(volumeSet, iBase_REGION);
222 : : #endif
223 [ + - ][ + - ]: 6 : me->commit_mesh(mit->second, COMPLETE_MESH);
224 : 6 : }
225 : :
226 : 6 : }
227 : :
228 : : // will also initialize NodeList, which comprises all <Vertex> data on the source sf
229 : 6 : void OneToOneSwept::BuildLateralBoundaryLayers(ModelEnt * me, std::vector<moab::EntityHandle> & layers)
230 : : {
231 : : // start with the source surface
232 : : //ModelEnt * sourceME = ModelEnt()
233 [ + - ]: 6 : iGeom * igeom_inst = me->igeom_instance();
234 : : iGeom::Error g_err;
235 : 6 : moab::ErrorCode rval = moab::MB_SUCCESS;
236 : : int index_id_src, index_id_tar;
237 : :
238 [ + - ]: 6 : g_err = igeom_inst->getIntData(sourceSurface, geom_id_tag, index_id_src);
239 [ + - ]: 6 : IBERRCHK(g_err, "Trouble get the int data for source surface.");
240 [ + - ]: 6 : g_err = igeom_inst->getIntData(targetSurface, geom_id_tag, index_id_tar);
241 [ + - ]: 6 : IBERRCHK(g_err, "Trouble get the int data for target surface.");
242 : :
243 : 6 : iBase_TagHandle taghandle = 0;
244 [ + - ][ + - ]: 6 : iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle);
[ + - ]
245 [ + + ]: 6 : if (m_err)
246 : : {
247 [ + - ][ + - ]: 3 : m_err = mk_core()->imesh_instance()->createTag("source",
248 [ + - ]: 3 : 1, iBase_INTEGER, taghandle);
249 : 3 : IBERRCHK(m_err,
250 [ + - ]: 3 : "Trouble creating the source tag handle in the mesh instance.");
251 : : }
252 : : // TODO: If the tag exists, we could verify that the tag is the proper one.
253 : : // Probably the tag should be put in a namespace, and perhaps it should
254 : : // be deleted when we are finished with it.
255 : :
256 [ + - ]: 6 : moab::Range quadsOnLateralSurfaces;
257 : 6 : ModelEnt * sME = NULL;
258 : :
259 [ + - ]: 12 : MEntVector surfs;
260 [ + - ]: 6 : me->get_adjacencies(2, surfs); // all surfaces adjacent to the volume
261 : :
262 [ + + ]: 42 : for (unsigned int i = 0; i < surfs.size(); i++)
263 : : {
264 : : int index_id_link;
265 [ + - ][ + - ]: 36 : g_err = igeom_inst->getIntData(surfs[i]->geom_handle(), geom_id_tag, index_id_link);
[ + - ]
266 [ + - ]: 36 : IBERRCHK(g_err, "Trouble get the int data for linking surface.");
267 [ + + ][ + + ]: 36 : if ((index_id_link != index_id_src) && (index_id_link != index_id_tar))
268 : : {
269 [ + - ][ + - ]: 24 : moab::EntityHandle surfSet = surfs[i]->mesh_handle();
270 : : // get all nodes from the surface, including the boundary nodes!
271 [ + - ]: 24 : rval = mb->get_entities_by_type(surfSet, moab::MBQUAD, quadsOnLateralSurfaces);
272 [ - + ][ # # ]: 24 : MBERRCHK(rval, mb);
[ # # ][ # # ]
273 : : }
274 [ + + ]: 36 : if (index_id_link == index_id_src)
275 [ + - ]: 6 : sME = surfs[i];
276 : : }
277 : : // the lateral layers will be build from these nodes
278 [ + - ]: 12 : moab::Range nodesOnLateralSurfaces;
279 [ + - ]: 6 : rval = mb->get_connectivity(quadsOnLateralSurfaces, nodesOnLateralSurfaces);
280 [ - + ][ # # ]: 6 : MBERRCHK(rval, mb);
[ # # ][ # # ]
281 : :
282 : : // construct NodeList, nodes on source surface
283 : : // list of node on boundary of source
284 [ + - ]: 12 : std::vector<moab::EntityHandle> bdyNodes;
285 [ + - ]: 12 : std::vector<int> group_sizes;
286 [ + - ]: 6 : sME-> boundary(0, bdyNodes, NULL, &group_sizes);// we do not need the senses, yet
287 : :
288 : 6 : numLoops = (int) group_sizes.size();
289 : 6 : sizeBLayer = (int) bdyNodes.size();
290 : :
291 [ + - ]: 6 : std::cout << "Execute one to one swept ....\n";
292 [ + - ][ + - ]: 6 : std::cout << "Nodes on boundary: " << bdyNodes.size() << "\n";
[ + - ]
293 [ + - ][ + - ]: 6 : std::cout << "Number of loops: " << group_sizes.size() << "\n";
[ + - ]
294 [ + - ][ + - ]: 6 : std::cout << "Nodes on lateral surfaces: " << nodesOnLateralSurfaces.size() << "\n";
[ + - ][ + - ]
295 [ + - ][ + - ]: 6 : std::cout << "Quads on lateral surfaces: " << quadsOnLateralSurfaces.size() << "\n";
[ + - ][ + - ]
296 : :
297 : 6 : numLayers = 0;
298 [ + - ]: 6 : if (!bdyNodes.empty())
299 [ + - ]: 6 : numLayers = nodesOnLateralSurfaces.size() / bdyNodes.size() - 1;
300 [ + - ][ - + ]: 6 : if (bdyNodes.size() * (numLayers + 1) != nodesOnLateralSurfaces.size())
301 : : {
302 [ # # ][ # # ]: 0 : std::cout << "Major problem: number of total nodes on boundary: " << nodesOnLateralSurfaces.size() << "\n";
[ # # ][ # # ]
303 [ # # ][ # # ]: 0 : std::cout << " number of nodes in one layer on the boundary:" << bdyNodes.size() << "\n";
[ # # ]
304 [ # # ][ # # ]: 0 : std::cout << " we expect bdyNodes.size()*(numLayers+1) == nodesOnLateralSurfaces.size(), but " << bdyNodes.size() * (numLayers + 1)
305 [ # # ][ # # ]: 0 : << " != " << nodesOnLateralSurfaces.size() << "\n";
[ # # ][ # # ]
306 : : }
307 : : // start from this list of boundary nodes (on the source)
308 : :
309 : : // get all nodes from the source surface
310 [ + - ]: 12 : moab::Range sourceNodes;
311 [ + - ]: 6 : moab::EntityHandle surfSet = sME->mesh_handle();
312 : : // get all nodes from the surface, including the boundary nodes!
313 [ + - ]: 12 : moab::Range quads;
314 [ + - ]: 6 : rval = mb->get_entities_by_type(surfSet, moab::MBQUAD, quads);
315 [ - + ][ # # ]: 6 : MBERRCHK(rval, mb);
[ # # ][ # # ]
316 : : // we do not really care about corners, only about boundary
317 [ + - ]: 6 : rval = mb->get_connectivity(quads, sourceNodes);
318 [ - + ][ # # ]: 6 : MBERRCHK(rval, mb);
[ # # ][ # # ]
319 : :
320 [ + - ][ + - ]: 6 : NodeList.resize(sourceNodes.size());
321 : 6 : int i = 0;
322 [ + - ][ + - ]: 1173 : for (moab::Range::iterator it = sourceNodes.begin(); it != sourceNodes.end(); it++, i++)
[ + - ][ + - ]
[ + + ]
323 : : {
324 [ + - ]: 1167 : moab::EntityHandle node = *it;
325 [ + - ][ + - ]: 1167 : NodeList[i].gVertexHandle = (iBase_EntityHandle) sourceNodes[i];
326 [ + - ]: 1167 : NodeList[i].index = i;
327 : :
328 [ + - ][ + - ]: 1167 : rval = mb->get_coords(&node, 1, &(NodeList[i].xyz[0]));
[ + - ]
329 [ - + ][ # # ]: 1167 : MBERRCHK(rval, mb);
[ # # ][ # # ]
330 [ + - ]: 1167 : NodeList[i].onBoundary = false;
331 : :
332 : : //set the int data to the entity
333 [ + - ][ + - ]: 1167 : m_err = mk_core()->imesh_instance()->setIntData(NodeList[i].gVertexHandle, taghandle, NodeList[i].index);
[ + - ][ + - ]
[ + - ]
334 [ + - ]: 1167 : IBERRCHK(m_err, "Trouble set the int value for nodes in the mesh instance.");
335 : : }
336 : : // now loop over the boundary and mark the boundary nodes as such
337 [ + + ]: 444 : for (unsigned int j = 0; j < bdyNodes.size(); j++)
338 : : {
339 [ + - ]: 438 : iBase_EntityHandle node = (iBase_EntityHandle) bdyNodes[j];
340 : : int index;
341 [ + - ][ + - ]: 438 : m_err = mk_core()->imesh_instance()->getIntData(node, taghandle, index);
[ + - ]
342 [ + - ]: 438 : IBERRCHK(m_err, "Trouble set the int value for nodes in the mesh instance.");
343 [ + - ]: 438 : NodeList[index].onBoundary = true;// some could be corners, but it is not that important
344 : : }
345 : :
346 : : // process faces on the source
347 [ + - ][ + - ]: 6 : FaceList.resize(quads.size());// quads is actually a range....
348 : :
349 [ + - ][ + + ]: 951 : for (unsigned int i = 0; i < quads.size(); i++)
350 : : {
351 [ + - ]: 945 : iBase_EntityHandle currentFace = (iBase_EntityHandle) quads[i];//
352 [ + - ]: 945 : FaceList[i].gFaceHandle = currentFace;
353 : : //FaceList[i].index = i;
354 : :
355 : : //get the nodes on the face elements
356 [ + - ]: 945 : std::vector<iBase_EntityHandle> Nodes;
357 [ + - ][ + - ]: 945 : m_err = mk_core()->imesh_instance()->getEntAdj(currentFace, iBase_VERTEX, Nodes);
[ + - ]
358 [ + - ]: 945 : IBERRCHK(m_err, "Trouble get the adjacent nodes from mesh face entity.");
359 : :
360 [ + - ][ + - ]: 945 : FaceList[i].connect.resize(Nodes.size());
361 [ + + ]: 4725 : for (unsigned int j = 0; j < Nodes.size(); j++)
362 : : {
363 : 3780 : int tmpIndex = -1;
364 [ + - ][ + - ]: 3780 : m_err = mk_core()->imesh_instance()->getIntData(Nodes[j], taghandle, tmpIndex);
[ + - ][ + - ]
365 [ + - ]: 3780 : IBERRCHK(m_err, "Trouble get the int data from node handle.");
366 : :
367 [ + - ][ + - ]: 3780 : FaceList[i].connect[j] = &NodeList[tmpIndex];// why not store tmp Index directly?
[ + - ]
368 : : }
369 [ + - ][ + - ]: 945 : m_err = mk_core()->imesh_instance()->setIntData(currentFace, taghandle, i);
[ + - ]
370 [ + - ]: 945 : IBERRCHK(m_err, "Trouble set the int data for quads on the source surface.");
371 : :
372 : : /* m_err = mk_core()->imesh_instance()->setIntData(currentFace, taghandle, i);
373 : : IBERRCHK(m_err, "Trouble set the int data for quad mesh on the source surface.");*/
374 : 945 : }
375 : :
376 [ + - ][ + - ]: 6 : std::copy(bdyNodes.begin(), bdyNodes.end(), std::back_inserter(layers));
377 : :
378 : : // mark all of the boundary nodes
379 [ + - ]: 6 : markedMoabEnts.insert(bdyNodes.begin(), bdyNodes.end());
380 : :
381 [ + + ]: 25 : for (int layer = 1; layer <= numLayers; layer++) // layer with index 0 is at the bottom/source
382 : : {
383 : 19 : int start_index = 0;//this is the start index in group
384 [ + + ]: 44 : for (unsigned int k = 0; k < group_sizes.size(); k++)
385 : : {
386 : : // first group starts at index 0, the rest are accumulated
387 [ + + ]: 25 : if (k > 0)
388 [ + - ]: 6 : start_index += group_sizes[k - 1];
389 [ + - ]: 25 : moab::EntityHandle node1 = layers[(layer - 1) * sizeBLayer + start_index];
390 [ + - ]: 25 : moab::EntityHandle node2 = layers[(layer - 1) * sizeBLayer + start_index + 1];
391 : : moab::EntityHandle node3, node4;
392 : : /*
393 : : * node3 -- node4 --> (next layer)
394 : : * | |
395 : : * node1 -- node2 ---> (guide layer)
396 : : */
397 [ + - ]: 25 : rval = NodeAbove(node1, node2, quadsOnLateralSurfaces, node3, node4);
398 [ - + ][ # # ]: 25 : MBERRCHK(rval, mb);
[ # # ][ # # ]
399 : : // check that node3 and node4 are not yet marked
400 [ + - ][ + - ]: 25 : if (markedMoabEnts.find(node3) != markedMoabEnts.end())
[ - + ]
401 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
402 [ + - ][ + - ]: 25 : if (markedMoabEnts.find(node4) != markedMoabEnts.end())
[ - + ]
403 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
404 [ + - ]: 25 : layers.push_back(node3);//start of a new layer
405 [ + - ]: 25 : layers.push_back(node4);//node above node 2
406 : : // mark node3 and node4
407 [ + - ]: 25 : markedMoabEnts.insert(node3);
408 [ + - ]: 25 : markedMoabEnts.insert(node4);
409 : : // we have at least 2 nodes in every loop, otherwise we are in deep trouble
410 : : //start
411 : 25 : int currentIndex = start_index + 2;
412 [ + - ][ + + ]: 1221 : while (currentIndex < start_index + group_sizes[k])
413 : : {
414 : 1196 : node1 = node2;
415 [ + - ]: 1196 : node2 = layers[(layer - 1) * sizeBLayer + currentIndex];// actually, the previous layer
416 : 1196 : node3 = node4;
417 [ + - ]: 1196 : rval = FourthNodeInQuad(node1, node2, node3, quadsOnLateralSurfaces, node4);
418 [ - + ][ # # ]: 1196 : MBERRCHK(rval, mb);
[ # # ][ # # ]
419 [ + - ][ + - ]: 1196 : if (markedMoabEnts.find(node4) != markedMoabEnts.end())
[ - + ]
420 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
421 [ + - ]: 1196 : layers.push_back(node4);
422 [ + - ]: 1196 : markedMoabEnts.insert(node4);
423 : 1196 : currentIndex++;
424 : : }
425 : : }// end group
426 : : }// end layer
427 [ + + ]: 31 : for (int lay = 0; lay <= numLayers; lay++)
428 : : {
429 [ + - ][ + - ]: 25 : std::cout << "layer : " << lay << "\n";
[ + - ]
430 [ + + ]: 1709 : for (int i = 0; i < sizeBLayer; i++)
431 : : {
432 [ + - ][ + - ]: 1684 : std::cout << " " << layers[lay * sizeBLayer + i];
[ + - ]
433 [ + + ]: 1684 : if (i % 10 == 9)
434 [ + - ]: 163 : std::cout << "\n";
435 : : }
436 [ + - ]: 25 : std::cout << "\n";
437 : 6 : }
438 : 6 : }
439 : :
440 : 25 : moab::ErrorCode OneToOneSwept::NodeAbove(moab::EntityHandle node1, moab::EntityHandle node2, moab::Range & latQuads,
441 : : moab::EntityHandle & node3, moab::EntityHandle & node4)
442 : : {
443 : : // look for node above node 1 in an unused quad
444 : 25 : moab::EntityHandle nds[2] = { node1, node2 };
445 : : // find all quads connected to these 2 nodes; find one in the range that is not used
446 [ + - ]: 25 : moab::Range adj_entities;
447 [ + - ]: 25 : moab::ErrorCode rval = mb->get_adjacencies(nds, 2, 2, false, adj_entities);
448 [ - + ][ # # ]: 25 : MBERRCHK(rval, mb);
[ # # ][ # # ]
449 : : // default is -1
450 : 25 : moab::EntityHandle nextQuad = 0;// null
451 [ + - ][ + - ]: 56 : for (unsigned int i = 0; i < adj_entities.size(); i++)
452 : : {
453 [ + - ][ + - ]: 152 : if ((markedMoabEnts.find(adj_entities[i]) == markedMoabEnts.end()) &&
[ + - ][ + + ]
[ + + ][ + - ]
[ + - ][ + - ]
[ + + # #
# # # # ]
454 [ + - ][ + - ]: 96 : (latQuads.find(adj_entities[i]) != latQuads.end()))
[ + - ][ + - ]
[ + + ][ + + ]
[ # # # # ]
455 : : {
456 : : // the quadrilateral is not marked but is on the lateral/linking surface
457 [ + - ]: 25 : nextQuad = adj_entities[i];
458 : 25 : break;
459 : : }
460 : : }
461 [ - + ]: 25 : if (0 == nextQuad)
462 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
463 : :
464 : : // decide on which side are nodes 1 and 2, then get the opposite side
465 : 25 : const moab::EntityHandle * conn4 = NULL;
466 : : int nnodes;
467 [ + - ]: 25 : rval = mb->get_connectivity(nextQuad, conn4, nnodes);
468 [ - + ][ # # ]: 25 : MBERRCHK(rval, mb);
[ # # ][ # # ]
469 : 25 : int index1 = -1;
470 [ + - ]: 37 : for (index1 = 0; index1 < 4; index1++)
471 [ + + ]: 37 : if (node1 == conn4[index1])
472 : 25 : break;
473 [ - + ]: 25 : if (4 == index1)
474 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
475 [ + + ]: 25 : if (node2 == conn4[(index1 + 1) % 4]) // quad is oriented node1, node2, node4, node3
476 : : {
477 : 19 : node3 = conn4[(index1 + 3) % 4];
478 : : }
479 [ + - ]: 6 : else if (node2 == conn4[(index1 +3) % 4]) // quad is oriented node2, node1, node3, node4
480 : : {
481 : 6 : node3 = conn4[(index1 + 1) % 4];
482 : : }
483 : : else
484 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);// something is really wrong
[ # # ]
485 : 25 : node4 = conn4[(index1 + 2) % 4];
486 : :
487 : : // mark the quad to be sure not to use it again
488 [ + - ]: 25 : markedMoabEnts.insert(nextQuad);
489 : :
490 : 25 : return moab::MB_SUCCESS;
491 : : }
492 : :
493 : 1196 : moab::ErrorCode OneToOneSwept::FourthNodeInQuad(moab::EntityHandle node1, moab::EntityHandle node2, moab::EntityHandle node3,
494 : : moab::Range & latQuads, moab::EntityHandle & node4)
495 : : {
496 : : // look for the fourth node
497 : 1196 : moab::EntityHandle nds[3] = { node1, node2, node3 };
498 : : // find all quads connected to these 3 nodes; find one in the range that is not used
499 [ + - ]: 1196 : moab::Range adj_entities;
500 [ + - ]: 1196 : moab::ErrorCode rval = mb->get_adjacencies(nds, 3, 2, false, adj_entities);
501 [ - + ][ # # ]: 1196 : MBERRCHK(rval, mb);
[ # # ][ # # ]
502 : :
503 : 1196 : moab::EntityHandle nextQuad = 0;// null
504 [ + - ][ + - ]: 1196 : for (unsigned int i = 0; i < adj_entities.size(); i++)
505 : : {
506 [ + - ][ + - ]: 3588 : if ((markedMoabEnts.find(adj_entities[i]) == markedMoabEnts.end()) &&
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - # #
# # # # ]
507 [ + - ][ + - ]: 2392 : (latQuads.find(adj_entities[i]) != latQuads.end()))
[ + - ][ + - ]
[ + - ][ + - ]
[ # # # # ]
508 : : {
509 : : // the quadrilateral is not marked but is on the lateral/linking surface
510 [ + - ]: 1196 : nextQuad = adj_entities[i];
511 : 1196 : break;
512 : : }
513 : : }
514 [ - + ]: 1196 : if (0 == nextQuad)
515 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
516 : :
517 : : // decide on which side are nodes 1 and 2, then get the opposite side
518 : 1196 : const moab::EntityHandle * conn4 = NULL;
519 : : int nnodes;
520 [ + - ]: 1196 : rval = mb->get_connectivity(nextQuad, conn4, nnodes);
521 [ - + ][ # # ]: 1196 : MBERRCHK(rval, mb);
[ # # ][ # # ]
522 : :
523 : 1196 : node4 = 0;
524 [ + - ]: 3100 : for (int index = 0; index < 4; index++)
525 : : {
526 : 3100 : moab::EntityHandle c4 = conn4[index];
527 [ + + ][ + + ]: 3100 : if (node1 != c4 && node2 != c4 && node3 != c4)
[ + + ]
528 : : {
529 : 1196 : node4 = c4;
530 : 1196 : break;
531 : : }
532 : : }
533 [ - + ]: 1196 : if (0 == node4)
534 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
535 : :
536 : : // mark the quad to be sure not to use it again
537 [ + - ]: 1196 : markedMoabEnts.insert(nextQuad);
538 : :
539 : 1196 : return moab::MB_SUCCESS;
540 : :
541 : : }
542 : : #ifdef HAVE_ARMADILLO
543 : : void OneToOneSwept::SurfMeshHarmonic(iBase_EntityHandle vol, vector<iBase_EntityHandle> &newNodehandle)
544 : : {
545 : :
546 : : SurfProHarmonicMap surf_pro(mk_core(), sourceSurface, targetSurface, vol);
547 : : surf_pro.match();
548 : : surf_pro.setMeshData(NodeList, TVertexList, FaceList);
549 : : surf_pro.projection();
550 : : surf_pro.getMeshData(TVertexList);
551 : :
552 : : iMesh::Error m_err;
553 : : for (unsigned int i = 0; i < TVertexList.size(); i++){
554 : : if (TVertexList[i].onBoundary)
555 : : continue;
556 : : m_err = mk_core()->imesh_instance()->createVtx(TVertexList[i].xyz[0], TVertexList[i].xyz[1], TVertexList[i].xyz[2], TVertexList[i].gVertexHandle);
557 : : IBERRCHK(m_err, "Trouble create a mesh nodes on the target surface!");
558 : : newNodehandle.push_back(TVertexList[i].gVertexHandle);
559 : : }
560 : : }
561 : : #endif
562 : :
563 : 0 : bool OneToOneSwept::isConcave()
564 : : {
565 : : //use the quad mesh on the source surface to determine whether it is concave or not
566 : 0 : iBase_TagHandle taghandle = 0;
567 [ # # ][ # # ]: 0 : iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle);
[ # # ]
568 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the tag handle in the mesh instance.");
569 [ # # ][ # # ]: 0 : for (std::vector<Vertex>::iterator it = NodeList.begin(); it != NodeList.end(); it++){
[ # # ]
570 [ # # ][ # # ]: 0 : if (!it->onBoundary)
571 : 0 : continue;
572 : : //get the adjacent quads around a vertex
573 : 0 : double angle = 0.0;
574 [ # # ]: 0 : std::vector<iBase_EntityHandle> adj_quads;
575 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getEntAdj(it->gVertexHandle, iBase_FACE, adj_quads);
[ # # ][ # # ]
576 [ # # ]: 0 : IBERRCHK(m_err, "Trouble get the adjacent quads around a vertex!");
577 [ # # ]: 0 : for (unsigned int i = 0; i < adj_quads.size(); i++){
578 : 0 : int face_index = -1;
579 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->getIntData(adj_quads[i], taghandle, face_index);
[ # # ][ # # ]
580 [ # # ]: 0 : if (m_err) continue;//this is a quad entity from the linking surface
581 : 0 : int node_index = 0, pre_index, next_index;
582 [ # # ]: 0 : for (int j = 0; j < 4; j++)
583 [ # # ][ # # ]: 0 : if (FaceList[face_index].connect[j]->index == it->index){
[ # # ][ # # ]
584 : 0 : node_index = j;
585 : 0 : break;
586 : : }
587 [ # # ][ # # ]: 0 : Vector3D v1(0.0), v2(0.0);
588 [ # # ][ # # ]: 0 : pre_index = (FaceList[face_index].getVertex((node_index+4-1)%4))->index;
589 [ # # ][ # # ]: 0 : next_index = (FaceList[face_index].getVertex((node_index+1)%4))->index;
590 [ # # ][ # # ]: 0 : v1 = NodeList[pre_index].xyz - it->xyz;
[ # # ][ # # ]
591 [ # # ][ # # ]: 0 : v2 = NodeList[next_index].xyz - it->xyz;
[ # # ][ # # ]
592 [ # # ][ # # ]: 0 : double dotproduct = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
[ # # ][ # # ]
[ # # ][ # # ]
593 [ # # ][ # # ]: 0 : angle += acos(dotproduct/sqrt((v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])*(v2[0]*v2[0]+v2[1]*v2[1]+v2[2]*v2[2])));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
594 : : }
595 [ # # ]: 0 : if (angle > M_PI){
596 [ # # ]: 0 : return true;
597 : : }
598 : 0 : }
599 : :
600 : 0 : return false;
601 : : }
602 : :
603 : :
604 : : //****************************************************************************//
605 : : // function : TargetSurfProjection
606 : : // Author : Shengyong Cai
607 : : // Date : Feb 15, 2011
608 : : // Description: map the mesh on the source surface to the target surface
609 : : //***************************************************************************//
610 : 6 : int OneToOneSwept::TargetSurfProjection(std::vector<moab::EntityHandle> & bLayers)
611 : : {
612 : : iMesh::Error m_err;
613 : : //iGeom::Error g_err;
614 : : iRel::Error r_err;
615 : : moab::ErrorCode rval;
616 : :
617 [ + - ]: 6 : MEntVector surfs;
618 [ + - ][ + - ]: 6 : mk_core()->get_entities_by_dimension(2, surfs);
619 : 6 : ModelEnt *target_surf = 0;
620 [ + - ]: 56 : for (unsigned int i = 0; i < surfs.size(); i++)
621 : : {
622 : :
623 [ + - ][ + - ]: 56 : if (surfs[i]->geom_handle() == targetSurface)
[ + + ]
624 : : {
625 [ + - ]: 6 : target_surf = surfs[i];
626 : 6 : break;
627 : : }
628 : : }
629 [ + - ]: 6 : int irelIndx = target_surf->iRelPairIndex();
630 [ + - ][ + - ]: 6 : iGeom * igeom_inst = mk_core()->igeom_instance(target_surf->iGeomIndex());
[ + - ]
631 : : // we could have got it from model tag, too
632 : : // something along this:
633 : : // err = igeom_instance(geom_index)->getData(*eit, iGeomModelTags[geom_index], &this_me);
634 : :
635 [ + - ]: 6 : TVertexList.resize(NodeList.size());
636 : : // make everything not on boundary, first; the true boundary flag will be set later
637 [ + + ]: 1173 : for (unsigned int i = 0; i < TVertexList.size(); i++)
638 : : {
639 [ + - ]: 1167 : TVertexList[i].onBoundary = false;
640 : : }
641 : :
642 : 6 : iBase_TagHandle taghandle_tar = 0;
643 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->getTagHandle("TargetMesh",
644 [ + - ]: 6 : taghandle_tar);
645 [ + + ]: 6 : if (m_err)
646 : : {
647 [ + - ][ + - ]: 3 : m_err = mk_core()->imesh_instance()->createTag("TargetMesh",
648 [ + - ]: 3 : 1, iBase_INTEGER, taghandle_tar);
649 : 3 : IBERRCHK(m_err,
650 [ + - ]: 3 : "Trouble create the target mesh tag handle in the mesh instance.");
651 : : }
652 : : // TODO: If the tag exists, we could verify that the tag is the proper one.
653 : : // Probably the tag should be put in a namespace, and perhaps it should
654 : : // be deleted when we are finished with it.
655 : :
656 : 6 : iBase_TagHandle taghandle = 0;
657 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle);
[ + - ]
658 [ + - ]: 6 : IBERRCHK(m_err, "Trouble get tag handle of the source surface.");
659 : :
660 : : // we know the first layer of nodes (on the source)
661 : : //
662 : : // the first 0 , 1, ..., sizeBlayer - 1 are on bottom, were bdyNodes before
663 : : //
664 [ + + ]: 444 : for (int i = 0; i < sizeBLayer; i++)
665 : : {
666 [ + - ]: 438 : iBase_EntityHandle sBNode = (iBase_EntityHandle) bLayers[i];
667 : : int index;
668 [ + - ][ + - ]: 438 : m_err = mk_core()->imesh_instance()->getIntData(sBNode, taghandle, index);
[ + - ]
669 [ + - ]: 438 : IBERRCHK(m_err, "Trouble get the int value for nodes in the mesh instance.");
670 [ + - ]: 438 : TVertexList[index].onBoundary = true;// some could be corners, but it is not that important
671 [ + - ]: 438 : moab::EntityHandle topNode = bLayers[numLayers * sizeBLayer + i];// node right above
672 : 438 : iBase_EntityHandle tBNode = (iBase_EntityHandle) topNode;// it is right above node i in loop
673 [ + - ]: 438 : TVertexList[index].gVertexHandle = tBNode;
674 : : // now get the coordinates of the tBNode (node on boundary of target)
675 [ + - ][ + - ]: 438 : rval = mb->get_coords(&topNode, 1, &(TVertexList[index].xyz[0]));
[ + - ]
676 [ - + ][ # # ]: 438 : MBERRCHK(rval, mb);
[ # # ][ # # ]
677 : :
678 : : }
679 : :
680 : : iBase_EntitySetHandle entityset; //this entityset is for storing the inner nodes on the target surface
681 [ + - ]: 12 : vector<iBase_EntityHandle> newNodehandle;
682 : :
683 [ + - ][ + - ]: 6 : r_err = mk_core()->irel_pair(irelIndx)->getEntSetRelation(targetSurface, 0, entityset);
[ + - ]
684 [ - + ]: 6 : if (r_err) //there is no entityset associated with targetSurface; this should be an error
685 : : {
686 [ # # ][ # # ]: 0 : m_err = mk_core()->imesh_instance()->createEntSet(1, entityset);
[ # # ]
687 [ # # ]: 0 : IBERRCHK(m_err, "Trouble create the entity set");
688 : : }
689 : :
690 : : // bool condition_harmonic = isConcave();
691 : 6 : bool is_exe_harmonic = false;
692 : : #ifdef HAVE_ARMADILLO
693 : : //target surface mesh projection based on harmonic mapping
694 : : if (condition_harmonic){
695 : : is_exe_harmonic = true;
696 : : SurfMeshHarmonic(volume, newNodehandle);
697 : : }
698 : : #endif
699 : :
700 [ + - ]: 6 : if (!is_exe_harmonic){
701 : : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
702 : : // Based on the transformation in the physical space, project the mesh nodes on the source surface to the target surface
703 : : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
704 [ + - ][ + - ]: 6 : Vector3D sPtsCenter(0.0), tPtsCenter(0.0);
705 [ + - ][ + - ]: 12 : std::vector<Vector3D> sBndNodes(sizeBLayer), tBndNodes(sizeBLayer);
706 : 6 : int num_pts = 0;
707 : : //calculate the barycenters for the source and target boundaries
708 [ + + ]: 1173 : for (unsigned int i = 0; i < NodeList.size(); i++)
709 : : {
710 [ + - ][ + + ]: 1167 : if (NodeList[i].onBoundary)
711 : : {
712 [ + - ][ + - ]: 438 : sPtsCenter += NodeList[i].xyz;
713 [ + - ][ + - ]: 438 : tPtsCenter += TVertexList[i].xyz;
714 [ + - ][ + - ]: 438 : sBndNodes[num_pts] = NodeList[i].xyz;
715 [ + - ][ + - ]: 438 : tBndNodes[num_pts] = TVertexList[i].xyz;
716 : 438 : num_pts++;
717 : : }
718 : : }
719 [ - + ]: 6 : if (sizeBLayer != num_pts)
720 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
721 [ + - ][ + - ]: 6 : sPtsCenter = sPtsCenter / num_pts;
722 [ + - ][ + - ]: 6 : tPtsCenter = tPtsCenter / num_pts;
723 : : //done with the barycenter calculation
724 : :
725 : : // compute transformation matrix from source boundary to target boundary
726 [ + - ]: 6 : Matrix3D transMatrix;
727 [ + - ]: 6 : computeTransformationFromSourceToTarget(sBndNodes, sPtsCenter, tBndNodes, tPtsCenter, transMatrix);
728 : :
729 : : //calculate the inner nodes on the target surface
730 [ + + ]: 1173 : for (unsigned int i = 0; i < NodeList.size(); i++)
731 : : {
732 [ + - ][ + + ]: 1167 : if (!NodeList[i].onBoundary)
733 : : {
734 [ + - ]: 729 : Vector3D xyz;
735 [ + - ][ + - ]: 729 : xyz = transMatrix * (NodeList[i].xyz - 2 * sPtsCenter + tPtsCenter)+sPtsCenter;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
736 : :
737 [ + - ][ + - ]: 729 : iGeom::Error g_err = igeom_inst->getEntClosestPt(targetSurface, xyz[0], xyz[1], xyz[2], TVertexList[i].xyz[0],
[ + - ][ + - ]
[ + - ]
738 [ + - ][ + - ]: 1458 : TVertexList[i].xyz[1], TVertexList[i].xyz[2]);
[ + - ][ + - ]
[ + - ]
739 [ + - ]: 729 : IBERRCHK(g_err, "Trouble get the closest point on the targets surface.");
740 : :
741 : : //create the node entities
742 [ + - ][ + - ]: 729 : iMesh::Error m_err = mk_core()->imesh_instance()->createVtx(TVertexList[i].xyz[0], TVertexList[i].xyz[1],
[ + - ][ + - ]
743 [ + - ][ + - ]: 1458 : TVertexList[i].xyz[2], TVertexList[i].gVertexHandle);
[ + - ][ + - ]
[ + - ][ + - ]
744 [ + - ]: 729 : IBERRCHK(m_err, "Trouble create the node entity.");
745 [ + - ]: 729 : TVertexList[i].onBoundary = false;
746 [ + - ][ + - ]: 729 : newNodehandle.push_back(TVertexList[i].gVertexHandle);
747 : : }
748 : 6 : }
749 : : }
750 : : //add the inner nodes to the entityset
751 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->addEntArrToSet(&newNodehandle[0], newNodehandle.size(), entityset);
[ + - ][ + - ]
752 [ + - ]: 6 : IBERRCHK(m_err, "Trouble add an array of nodes to the entityset.");
753 : :
754 : : //until now, all the nodes have been generated on the target surface
755 : :
756 : : // we need to decide the orientation with respect to the faces on the source (which are oriented correctly)
757 : : // look at the first 2 nodes in the boundary, in target
758 : 6 : int sense_out = 1;
759 [ + - ]: 12 : std::vector<moab::EntityHandle> bdyTNodes;
760 [ + - ]: 12 : std::vector<int> group_sizes;
761 [ + - ]: 6 : target_surf-> boundary(0, bdyTNodes, NULL, &group_sizes);// we do not need the senses, yet. Or do we?
762 : : // find the first node and second node in the bLayers[] list corresponding to the top
763 [ + - ][ + - ]: 6 : moab::EntityHandle node1 = bdyTNodes[1], node2 = bdyTNodes[2];
764 : 6 : int index = -1;
765 [ + - ]: 154 : for (index = 0; index < sizeBLayer; index++)
766 : : {
767 [ + - ][ + + ]: 154 : if (bLayers[numLayers * sizeBLayer + index] == node1)
768 : 6 : break;
769 : : }
770 [ - + ]: 6 : if (index == sizeBLayer)
771 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);
[ # # ]
772 : :
773 : 6 : int prevIndex = (index + sizeBLayer -1) % sizeBLayer;
774 : 6 : int nextIndex = (index + 1) % sizeBLayer;
775 [ + - ][ - + ]: 6 : if (bLayers[numLayers * sizeBLayer + prevIndex] == node2)
776 : 0 : sense_out = -1;
777 [ + - ][ + - ]: 6 : else if (bLayers[numLayers * sizeBLayer + nextIndex] == node2)
778 : 6 : sense_out = 1;
779 : : else
780 [ # # ][ # # ]: 0 : MBERRCHK(moab::MB_FAILURE, mb);// serious error in the logic, can't decide orientation
[ # # ]
781 : :
782 : :
783 : : //create the quadrilateral elements on the target surface
784 [ + - ]: 12 : vector<iBase_EntityHandle> mFaceHandle(FaceList.size());
785 [ + - ]: 12 : vector<iBase_EntityHandle> connect(FaceList.size() * 4);
786 [ + + ]: 951 : for (unsigned int i = 0; i < FaceList.size(); i++)
787 : : {
788 [ - + ]: 945 : if (sense_out < 0)
789 : : {
790 [ # # ][ # # ]: 0 : connect[4 * i + 0] = TVertexList[(FaceList[i].getVertex(0))->index].gVertexHandle;
[ # # ][ # # ]
791 [ # # ][ # # ]: 0 : connect[4 * i + 1] = TVertexList[(FaceList[i].getVertex(3))->index].gVertexHandle;
[ # # ][ # # ]
792 [ # # ][ # # ]: 0 : connect[4 * i + 2] = TVertexList[(FaceList[i].getVertex(2))->index].gVertexHandle;
[ # # ][ # # ]
793 [ # # ][ # # ]: 0 : connect[4 * i + 3] = TVertexList[(FaceList[i].getVertex(1))->index].gVertexHandle;
[ # # ][ # # ]
794 : : }
795 : : else
796 : : {
797 [ + - ][ + - ]: 945 : connect[4 * i + 0] = TVertexList[(FaceList[i].getVertex(0))->index].gVertexHandle;
[ + - ][ + - ]
798 [ + - ][ + - ]: 945 : connect[4 * i + 1] = TVertexList[(FaceList[i].getVertex(1))->index].gVertexHandle;
[ + - ][ + - ]
799 [ + - ][ + - ]: 945 : connect[4 * i + 2] = TVertexList[(FaceList[i].getVertex(2))->index].gVertexHandle;
[ + - ][ + - ]
800 [ + - ][ + - ]: 945 : connect[4 * i + 3] = TVertexList[(FaceList[i].getVertex(3))->index].gVertexHandle;
[ + - ][ + - ]
801 : : }
802 : :
803 : : }
804 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->createEntArr(iMesh_QUADRILATERAL, &connect[0], connect.size(), &mFaceHandle[0]);
[ + - ][ + - ]
[ + - ]
805 [ + - ]: 6 : IBERRCHK(m_err, "Trouble create the quadrilateral mesh.");
806 : :
807 : : //add the inner face elements to the entityset
808 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->addEntArrToSet(&mFaceHandle[0], FaceList.size(), entityset);
[ + - ][ + - ]
809 [ + - ]: 6 : IBERRCHK(m_err, "Trouble add an array of quadrilateral entities to the entity set.");
810 : :
811 [ + - ][ + - ]: 6 : mk_core()->save_mesh("BeforeHex.vtk");
812 : : #ifdef HAVE_MESQUITE
813 : :
814 [ + - ]: 6 : SurfMeshOptimization();
815 : :
816 : : #endif
817 : :
818 : 6 : return 1;
819 : : }
820 : : // input: list of nodes on source, boundary center, list of nodes on target, target center
821 : : // output: 3x3 matrix A such that
822 : : // target= A * ( source - 2*sc + tc) + sc
823 : 32 : void OneToOneSwept::computeTransformationFromSourceToTarget(std::vector<Vector3D> & sNodes, Vector3D & sc,
824 : : std::vector<Vector3D> & tNodes, Vector3D & tc, Matrix3D & transMatrix)
825 : : {
826 [ + - ]: 32 : Matrix3D tmpMatrix(0.0);
827 [ + - ]: 32 : Matrix3D bMatrix(0.0);
828 : : //transform the coordinates
829 : 32 : int size = (int)sNodes.size();
830 [ - + ]: 32 : assert(sNodes.size() == tNodes.size());
831 [ + + ]: 2086 : for (int i = 0; i < size; i++)
832 : : {
833 : : //transform the boundary nodes
834 [ + - ][ + - ]: 2054 : sNodes[i] = sNodes[i] - 2 * sc + tc;
[ + - ][ + - ]
[ + - ][ + - ]
835 [ + - ][ + - ]: 2054 : tNodes[i] = tNodes[i] - sc;
[ + - ][ + - ]
836 : : }
837 : :
838 : : //calculate the transformation matrix: transform the nodes on the source surface to the target surface in the physical space
839 [ + + ]: 2086 : for (int i = 0; i < size; i++)
840 : : {
841 [ + - ][ + - ]: 2054 : tmpMatrix = tmpMatrix + sNodes[i]*transpose(sNodes[i]);
[ + - ][ + - ]
[ + - ]
842 [ + - ][ + - ]: 2054 : bMatrix = bMatrix + tNodes[i]*transpose(sNodes[i]);
[ + - ][ + - ]
[ + - ]
843 : : }
844 : :
845 : : //first determine the affine mapping matrix is singular or not
846 [ + - ]: 32 : double detValue = det(tmpMatrix);
847 : : (void) detValue;
848 [ + - ][ - + ]: 32 : assert(pow(detValue, 2)>1.0e-20);
849 : :
850 : : //solve the affine mapping matrix, make use of inverse matrix to get affine mapping matrix
851 [ + - ]: 32 : Matrix3D InvMatrix = inverse(tmpMatrix);
852 [ + - ]: 32 : transMatrix = bMatrix*InvMatrix;
853 : :
854 : 32 : }
855 : : // use similar code to TargetSurfProjection, but do not project on surface...
856 : 6 : int OneToOneSwept::ProjectInteriorLayers(std::vector<moab::EntityHandle> & boundLayers, vector<vector<Vertex> > & linkVertexList)
857 : : {
858 : : iMesh::Error m_err;
859 [ + - ][ + - ]: 6 : Vector3D sPtsCenter(0.), tPtsCenter(0.);
860 [ + - ]: 6 : std::vector<Vector3D> PtsCenter(numLayers - 1);
861 [ + - ][ + - ]: 12 : std::vector<Vector3D> sBoundaryNodes(0), tBoundaryNodes(0);
862 : :
863 : 6 : iBase_TagHandle taghandle = 0;
864 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->getTagHandle("source", taghandle);
[ + - ]
865 [ + - ]: 6 : IBERRCHK(m_err, "Trouble getting the source tag handle");
866 : :
867 [ + - ]: 12 : std::vector<Vector3D> latCoords;
868 [ + - ]: 6 : latCoords.resize(boundLayers.size() - 2 * sizeBLayer);// to store the coordinates of the lateral points
869 : : // (exclude the source and target, they are already in NodeList and TVertexList)...
870 : :
871 [ + - ][ + - ]: 6 : moab::ErrorCode rval = mb->get_coords(&(boundLayers[sizeBLayer]), latCoords.size(), &(latCoords[0][0])); //these are the nodes for lateral sides
[ + - ][ + - ]
872 [ - + ][ # # ]: 6 : MBERRCHK(rval, mb);
[ # # ][ # # ]
873 : :
874 : : //calculate the center coordinates
875 [ + + ]: 19 : for (int i = 0; i < numLayers - 1; i++)
876 [ + - ][ + - ]: 13 : PtsCenter[i].set(0.);
877 : :
878 : 6 : int numPts = sizeBLayer;// a lot of repetition ; do we really need
879 [ + - ]: 6 : sBoundaryNodes.resize(numPts);// another useless copy
880 [ + - ]: 6 : tBoundaryNodes.resize(numPts);// another useless copy
881 [ + + ]: 444 : for (int j = 0; j < sizeBLayer; j++)
882 : : {
883 [ + - ]: 438 : iBase_EntityHandle sBNode = (iBase_EntityHandle) boundLayers[j];
884 : : int i;
885 : : // this is the index in NodeList...
886 [ + - ][ + - ]: 438 : m_err = mk_core()->imesh_instance()->getIntData(sBNode, taghandle, i);
[ + - ]
887 [ + - ]: 438 : IBERRCHK(m_err, "Trouble get the int value for nodes in the mesh instance.");
888 : :
889 [ + - ][ + - ]: 438 : sPtsCenter += NodeList[i].xyz;
890 [ + - ][ + - ]: 438 : tPtsCenter += TVertexList[i].xyz;
891 : :
892 [ + - ][ + - ]: 438 : sBoundaryNodes[j] = NodeList[i].xyz;
893 [ + - ][ + - ]: 438 : tBoundaryNodes[j] = TVertexList[i].xyz;
894 : :
895 [ + + ]: 1246 : for (int lay = 0; lay < numLayers - 1; lay++)
896 : : {
897 : : // interior layers
898 : 808 : int indexNodeOnSide = lay * sizeBLayer + j;// all p3d will be above index i
899 [ + - ]: 808 : Vector3D & p3d = latCoords[indexNodeOnSide];
900 [ + - ][ + - ]: 808 : PtsCenter[lay] += p3d;
901 [ + - ][ + - ]: 808 : linkVertexList[lay][i].xyz = p3d;
902 [ + - ][ + - ]: 808 : linkVertexList[lay][i].gVertexHandle = (iBase_EntityHandle) boundLayers[indexNodeOnSide + sizeBLayer];
[ + - ]
903 : : }
904 : : }
905 : :
906 [ + - ][ + - ]: 6 : sPtsCenter = sPtsCenter / numPts;
907 [ + - ][ + - ]: 6 : tPtsCenter = tPtsCenter / numPts;
908 : :
909 : : //calculate the center coordinates for the ith layer
910 [ + + ]: 19 : for (int i = 0; i < numLayers - 1; i++)
911 [ + - ][ + - ]: 13 : PtsCenter[i] = PtsCenter[i] / numPts;
[ + - ][ + - ]
912 : :
913 [ + - ]: 12 : std::vector<Vector3D> sNodes(numPts); //boundary nodes on the source surface
914 [ + - ]: 12 : std::vector<Vector3D> tNodes(numPts); //boundary nodes on the target surface
915 [ + - ][ + - ]: 12 : std::vector<Vector3D> isBNodes(numPts), itBNodes(numPts); //boundary nodes on the different layers
916 : :
917 : : //loop over different layers
918 [ + + ]: 19 : for (int i = 0; i < numLayers - 1; i++)
919 : : {
920 [ + + ]: 821 : for (int j = 0; j < numPts; j++)
921 : : {
922 [ + - ][ + - ]: 808 : sNodes[j] = sBoundaryNodes[j];
923 [ + - ][ + - ]: 808 : tNodes[j] = tBoundaryNodes[j];
924 : : //coordinates on the different layers
925 [ + - ][ + - ]: 808 : isBNodes[j] = itBNodes[j] = latCoords[i*sizeBLayer + j];
[ + - ]
926 : : }
927 : :
928 [ + - ][ + - ]: 13 : Matrix3D sA, tA;
929 : : // from source to layer i
930 [ + - ][ + - ]: 13 : computeTransformationFromSourceToTarget(sNodes, sPtsCenter, isBNodes, PtsCenter[i], sA);
931 : : // from target to layer i
932 [ + - ][ + - ]: 13 : computeTransformationFromSourceToTarget(tNodes, tPtsCenter, itBNodes, PtsCenter[i], tA);
933 : : //calculate the inner nodes for different layers
934 : 13 : double s = (i + 1) / double(numLayers); // interpolation factor
935 [ + + ]: 2258 : for (unsigned int j = 0; j < NodeList.size(); j++)
936 : : {
937 [ + - ][ + + ]: 2245 : if (!(NodeList[j].onBoundary))
938 : : {
939 [ + - ][ + - ]: 1437 : Vector3D spts, tpts, pts;
[ + - ]
940 : : iBase_EntityHandle nodeHandle;
941 [ + - ][ + - ]: 1437 : spts = sA * (NodeList[j].xyz - 2 * sPtsCenter + PtsCenter[i])+sPtsCenter;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
942 [ + - ][ + - ]: 1437 : tpts = tA * (TVertexList[j].xyz - 2 * tPtsCenter + PtsCenter[i])+tPtsCenter;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
943 : : // interpolate the 2 results
944 [ + - ][ + - ]: 1437 : pts = (1-s) * spts + s * tpts;
[ + - ][ + - ]
945 : :
946 [ + - ][ + - ]: 1437 : linkVertexList[i][j].xyz = pts;
947 : :
948 [ + - ][ + - ]: 1437 : m_err = mk_core()->imesh_instance()->createVtx(pts[0], pts[1], pts[2], nodeHandle);
[ + - ][ + - ]
[ + - ][ + - ]
949 [ + - ]: 1437 : IBERRCHK(m_err, "Trouble create the vertex entity.");
950 [ + - ][ + - ]: 1437 : linkVertexList[i][j].gVertexHandle = nodeHandle;
951 [ + - ][ + - ]: 1437 : m_err = mk_core()->imesh_instance()->addEntToSet(nodeHandle, volumeSet);
[ + - ]
952 [ + - ]: 1437 : IBERRCHK(m_err, "Trouble add the node handle to the entity set.");
953 : : }
954 : : }
955 : : }
956 : 6 : return 1;
957 : : }
958 : :
959 : : //****************************************************************************//
960 : : // function : CreateElements
961 : : // Author : Shengyong Cai
962 : : // Date : Feb 16, 2011
963 : : // Description: create hexahedral elements by connecting 8 nodes which have
964 : : // already been created by previous functions
965 : : //***************************************************************************//
966 : 6 : int OneToOneSwept::CreateElements(vector<vector<Vertex> > &linkVertexList)
967 : : {
968 : : //create the quadrilateral elements on the different layers
969 : : //it is not necessary to create the quadrilateral elements on the different layers. Hex element can be created based on the eight nodes
970 : : iMesh::Error m_err;
971 [ + - ]: 6 : vector<iBase_EntityHandle> mVolumeHandle(FaceList.size());
972 : : // first decide orientation, based on the FaceList orientation, and first node above
973 : : // take one face on source, and first node above in layer 1
974 [ + - ][ + - ]: 6 : Vertex * v1 = FaceList[0].connect[0];
975 [ + - ][ + - ]: 6 : Vertex * v2 = FaceList[0].connect[1];
976 [ + - ][ + - ]: 6 : Vertex * v3 = FaceList[0].connect[2];
977 : : // vertex 4 is on layer 1, above vertex 1
978 : 6 : int ind1 = v1->index;
979 : : Vertex * v4;
980 [ + - ]: 6 : if (numLayers > 1)
981 [ + - ][ + - ]: 6 : v4 = &linkVertexList[0][ind1];
982 : : else
983 [ # # ]: 0 : v4 = &TVertexList[ind1];
984 : :
985 [ + - ][ + - ]: 6 : Vector3D normal1 = (v2->xyz-v1->xyz)*(v3->xyz-v1->xyz);
[ + - ][ + - ]
986 [ + - ][ + - ]: 6 : double vol6= normal1 % (v4->xyz-v1->xyz);
987 [ + - ]: 6 : std::cout << "Orientation is ";
988 : 6 : int skip = 0;
989 [ + + ]: 6 : if (vol6 < 0)
990 : : {
991 [ + - ]: 5 : std::cout <<"negative\n";
992 : 5 : skip=4;
993 : : }
994 : : else
995 [ + - ]: 1 : std::cout <<"positive\n";
996 : :
997 [ + - ]: 12 : vector<iBase_EntityHandle> connect(8);
998 [ + + ]: 25 : for (int m = 0; m < numLayers; m++)
999 : : {
1000 [ + + ]: 2795 : for (unsigned int i = 0; i < FaceList.size(); i++)
1001 : : {
1002 [ + + ]: 13880 : for (int k = 0; k < 4; k++)
1003 : : {
1004 [ + + ]: 11104 : if (m == 0) // the first layer is the source layer
1005 [ + - ]: 3780 : connect[k + skip] = NodeList
1006 [ + - ][ + - ]: 3780 : [(FaceList[i].connect[k])->index].gVertexHandle;
[ + - ]
1007 : : else // the first layer is an intermediate layer
1008 [ + - ][ + - ]: 14648 : connect[k + skip] = linkVertexList[m - 1]
1009 [ + - ][ + - ]: 14648 : [(FaceList[i].connect[k])->index].gVertexHandle;
[ + - ]
1010 [ + + ]: 11104 : if (m == (numLayers - 1)) // the second layer is the target layer
1011 [ + - ]: 3780 : connect[(k + skip + 4) % 8] = TVertexList
1012 [ + - ][ + - ]: 3780 : [(FaceList[i].connect[k])->index].gVertexHandle;
[ + - ]
1013 : : else // the second layer is an intermediate layer
1014 [ + - ][ + - ]: 14648 : connect[(k + skip + 4) % 8] = linkVertexList[m]
1015 [ + - ][ + - ]: 14648 : [(FaceList[i].connect[k])->index].gVertexHandle;
[ + - ]
1016 : : }
1017 : : m_err = mk_core()->imesh_instance()->createEnt(iMesh_HEXAHEDRON,
1018 [ + - ][ + - ]: 2776 : &connect[0], 8, mVolumeHandle[i]);
[ + - ][ + - ]
[ + - ]
1019 [ + - ]: 2776 : IBERRCHK(m_err, "Trouble creating the hexahedral element.");
1020 : : }
1021 [ + - ][ + - ]: 19 : m_err = mk_core()->imesh_instance()->addEntArrToSet(&mVolumeHandle[0],
[ + - ]
1022 [ + - ]: 38 : FaceList.size(), volumeSet);
1023 : 19 : IBERRCHK(m_err,
1024 [ + - ]: 19 : "Trouble adding an array of hexahedral elements to the entity set.");
1025 : : }
1026 : 6 : return 1;
1027 : : }
1028 : :
1029 : : #ifdef HAVE_MESQUITE
1030 : : //target surface mesh smoothing by Mesquite
1031 : 6 : void OneToOneSwept::SurfMeshOptimization()
1032 : : {
1033 : 6 : MEntSelection::iterator mit = mentSelection.begin();
1034 [ + - ]: 6 : ModelEnt *me = mit -> first;
1035 [ + - ]: 6 : int irelPairIndex = me->iRelPairIndex();
1036 [ + - ][ + - ]: 6 : iGeom * igeom_inst = mk_core()->igeom_instance(me->iGeomIndex());
[ + - ]
1037 : : //create a tag to attach the coordinates to nodes
1038 : 6 : iBase_TagHandle mapped_tag = 0;
1039 [ + - ][ + - ]: 6 : iMesh::Error m_err = mk_core()->imesh_instance()->getTagHandle("MsqAltCoords", mapped_tag);
[ + - ]
1040 [ + - ]: 6 : if (m_err)
1041 : : {
1042 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->createTag("MsqAltCoords", 3, iBase_DOUBLE, mapped_tag);
[ + - ]
1043 [ + - ]: 6 : IBERRCHK(m_err, "Trouble create a tag.");
1044 : :
1045 : : }
1046 : : //attach the coordinates to the nodes
1047 : 6 : double tag_data[3*TVertexList.size()];
1048 [ + - ]: 6 : std::vector<iBase_EntityHandle> vertexHandle(TVertexList.size());
1049 [ + + ]: 1173 : for (unsigned int i = 0; i < NodeList.size();i++)
1050 : : {
1051 [ + - ][ + - ]: 1167 : tag_data[3*i] = NodeList[i].xyz[0];
1052 [ + - ][ + - ]: 1167 : tag_data[3*i+1] = NodeList[i].xyz[1];
1053 [ + - ][ + - ]: 1167 : tag_data[3*i+2] = NodeList[i].xyz[2];
1054 [ + - ][ + - ]: 1167 : vertexHandle[i] = TVertexList[i].gVertexHandle;
1055 : : }
1056 [ + - ][ + - ]: 6 : m_err = mk_core()->imesh_instance()->setDblArrData(&vertexHandle[0], NodeList.size(), mapped_tag, &tag_data[0]);
[ + - ][ + - ]
1057 [ + - ]: 6 : IBERRCHK(m_err, "Trouble set an array of int data to nodes.");
1058 : :
1059 : : //get the mesh entityset for target surface
1060 : : iBase_EntitySetHandle surfSets;
1061 [ + - ][ + - ]: 6 : iRel::Error r_err = mk_core()->irel_pair(irelPairIndex)->getEntSetRelation(targetSurface, 0, surfSets);
[ + - ]
1062 [ + - ]: 6 : IBERRCHK(r_err, "Trouble get the mesh entity set for the target surface.");
1063 : : //call the MeshImprove class to smooth the target surface mesh by using Mesquite
1064 [ + - ][ + - ]: 12 : MeshImprove meshopt(mk_core(), false, true, true, true, igeom_inst);
1065 [ + - ]: 12 : meshopt.SurfMeshImprove(targetSurface, surfSets, iBase_FACE);
1066 : 6 : }
1067 : : #endif
1068 [ + - ][ + - ]: 156 : }
1069 : :
|
