Branch data Line data Source code
1 : : #include "meshkit/EdgeMesher.hpp"
2 : : #include "meshkit/Matrix.hpp"
3 : : #include "meshkit/MKCore.hpp"
4 : : #include "meshkit/ModelEnt.hpp"
5 : : #include "meshkit/SizingFunction.hpp"
6 : : #include "meshkit/RegisterMeshOp.hpp"
7 : : #include "moab/ReadUtilIface.hpp"
8 : : #include <vector>
9 : : #include <math.h>
10 : :
11 : : namespace MeshKit
12 : : {
13 : :
14 : : //---------------------------------------------------------------------------//
15 : : //Entity Type initilization for edge meshing
16 : : moab::EntityType EdgeMesher_tps[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBMAXTYPE};
17 : 40 : const moab::EntityType* EdgeMesher::output_types()
18 : 40 : { return EdgeMesher_tps; }
19 : :
20 : : //---------------------------------------------------------------------------//
21 : : // Construction Function for Edge Mesher
22 : 60 : EdgeMesher::EdgeMesher(MKCore *mk_core, const MEntVector &me_vec)
23 : 60 : : MeshScheme(mk_core, me_vec), schemeType(EQUAL), ratio(1.2)
24 : : {
25 : :
26 : 60 : }
27 : : #if 0
28 : : //---------------------------------------------IBERRCHK(g_err, "Trouble get the adjacent geometric nodes on a surface.");------------------------------//
29 : : // measure function: compute the distance between the parametric coordinate
30 : : // ustart and the parametric coordinate uend
31 : : double EdgeMesher::measure(iGeom::EntityHandle ent, double ustart, double uend) const
32 : : {
33 : : double umin, umax;
34 : :
35 : : //get the minimal and maximal parametrical coordinates of the edge
36 : : iGeom::Error err = mk_core()->igeom_instance()->getEntURange(ent, umin, umax);
37 : : IBERRCHK(err, "Trouble getting parameter range for edge.");
38 : :
39 : : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to same parameter umax and umin.");
40 : :
41 : : //compute the distance for edges
42 : : double measure;
43 : : err = mk_core()->igeom_instance()->measure(&ent, 1, &measure);
44 : :
45 : : IBERRCHK(err, "Trouble getting edge measure.");
46 : :
47 : : return measure * (uend - ustart) / (umax - umin);
48 : : }
49 : : #endif
50 : : //---------------------------------------------------------------------------//
51 : : // setup function: set up the number of intervals for edge meshing through the
52 : : // sizing function
53 : 61 : void EdgeMesher::setup_this()
54 : : {
55 : : //compute the number of intervals for the associated ModelEnts, from the size set on them
56 : : //the sizing function they point to, or a default sizing function
57 [ + - ][ + - ]: 352 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
58 : : {
59 [ + - ]: 291 : ModelEnt *me = mit->first;
60 : :
61 : : //first check to see whether entity is meshed
62 [ + - ][ + + ]: 291 : if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0)
[ + - ][ + + ]
[ + + ]
63 : 32 : continue;
64 : :
65 [ + - ][ + - ]: 259 : SizingFunction *sf = mk_core()->sizing_function(me->sizing_function_index());
[ + - ]
66 [ + + ][ + - ]: 263 : if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT &&
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ]
67 [ + - ][ + - ]: 4 : mk_core()->sizing_function(0))
68 : : {
69 [ + - ][ + - ]: 4 : sf = mk_core()->sizing_function(0);
70 [ + - ]: 4 : me->sizing_function_index(0);
71 : : }
72 : :
73 [ - + ][ # # ]: 259 : if (!sf && me -> mesh_intervals() < 0 && me -> interval_firmness() == DEFAULT)
[ # # ][ # # ]
[ # # ][ - + ]
74 : : {
75 : : //no sizing set, just assume default #intervals as 4
76 [ # # ]: 0 : me->mesh_intervals(4);
77 [ # # ]: 0 : me->interval_firmness(DEFAULT);
78 : : }
79 : : else
80 : : {
81 : : //check # intervals first, then size, and just choose for now
82 [ + - ][ + + ]: 259 : if (sf->intervals() > 0)
83 : : {
84 [ + - ][ - + ]: 32 : if (me->constrain_even() && sf->intervals()%2)
[ # # ][ # # ]
[ - + ]
85 [ # # ][ # # ]: 0 : me -> mesh_intervals(sf->intervals()+1);
86 : : else
87 [ + - ][ + - ]: 32 : me -> mesh_intervals(sf->intervals());
88 [ + - ]: 32 : me -> interval_firmness(HARD);
89 : : }
90 [ + - ][ + - ]: 227 : else if (sf->size()>0)
91 : : {
92 [ + - ][ + - ]: 227 : int intervals = me->measure()/sf->size();
93 [ - + ]: 227 : if (!intervals) intervals++;
94 [ + - ][ + + ]: 227 : if (me->constrain_even() && intervals%2) intervals++;
[ + + ][ + + ]
95 [ + - ]: 227 : me->mesh_intervals(intervals);
96 [ + - ]: 227 : me->interval_firmness(SOFT);
97 : : }
98 : : else
99 [ # # ]: 0 : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION, "Sizing function for edge had neither positive size nor positive intervals.");
100 : : }
101 : : }
102 : :
103 : : // now call setup_boundary to treat vertices
104 : : // Wrong!!!!!!!!!
105 : : // setup_boundary();
106 : : /* this is not enough to ensure that vertex mesher will be called before
107 : : "this" edge mesher
108 : : the case that fell through the cracks was if the end nodes were already setup
109 : : then the this_op[0] would not be retrieved, and not inserted "before" the edge mesher MeshOp
110 : : */
111 : 61 : int dim=0;
112 [ + - ][ + - ]: 61 : MeshOp * vm = (MeshOp*) mk_core()->construct_meshop(dim);
113 : :
114 [ + - ][ + - ]: 352 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
115 : : {
116 [ + - ]: 291 : ModelEnt *me = mit->first;
117 [ + - ]: 291 : MEntVector children;
118 [ + - ]: 291 : me->get_adjacencies(0, children);
119 [ + + ]: 857 : for (unsigned int i=0; i<children.size(); i++)
120 [ + - ][ + - ]: 566 : if (children[i]->is_meshops_list_empty())
[ + + ]
121 : : {
122 [ + - ][ + - ]: 222 : vm->add_modelent(children[i]);
123 : : }
124 : 291 : }
125 : : // in any case, make sure that the vertex mesher is inserted before this edge mesher
126 : 61 : mk_core()->insert_node(vm, this, mk_core()->root_node());
127 : :
128 : :
129 : 61 : }
130 : :
131 : : //---------------------------------------------------------------------------//
132 : : // execute function: Generate the mesh for edges
133 : 61 : void EdgeMesher::execute_this()
134 : : {
135 [ + - ]: 61 : std::vector<double> coords;
136 [ + - ]: 122 : std::vector<moab::EntityHandle> nodes;
137 : :
138 [ + - ][ + - ]: 352 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
139 : : {
140 [ + - ]: 291 : ModelEnt *me = mit -> first;
141 [ + - ][ + + ]: 291 : if (me->get_meshed_state() >= COMPLETE_MESH)
142 : 4 : continue;
143 : : //resize the coords based on the interval setting
144 [ + - ]: 287 : int num_edges = me->mesh_intervals();
145 [ + - ]: 287 : coords.resize(3*(num_edges+1));
146 : 287 : nodes.clear();
147 [ + - ]: 287 : nodes.reserve(num_edges + 1);
148 : :
149 : : //get bounding mesh entities, use 1st 2 entries of nodes list temporarily
150 : : //pick up the boundary end nodes
151 [ + - ]: 287 : me->boundary(0, nodes);
152 : :
153 : 287 : bool periodic = (nodes.size() == 1);
154 : :
155 : : //get coords in list, then move one tuple to the last position
156 [ + - ][ + - ]: 287 : moab::ErrorCode rval = mk_core()->moab_instance()->get_coords(&nodes[0], nodes.size(), &coords[0]);
[ + - ][ + - ]
[ + - ]
157 [ - + ][ # # ]: 287 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
158 : :
159 : : //move the second node to the endmost position in the node list
160 : : // if periodic, the last node coordinates are also at index 0 in coords array
161 : : // there is only one node, coords[3+i] are not even initialized
162 [ + + ]: 287 : int index2 = (periodic) ? 0 : 3;
163 [ + + ]: 1148 : for (int i = 0; i < 3; i++)
164 [ + - ][ + - ]: 861 : coords[3*num_edges+i] = coords[index2+i];
165 : :
166 : 287 : EdgeSchemeType scheme = schemeType;
167 [ + - ][ + - ]: 287 : SizingFunction *sz = mk_core()->sizing_function(me->sizing_function_index());
[ + - ]
168 [ + - ][ + + ]: 287 : if (sz->variable())
169 : 19 : scheme = VARIABLE;
170 : :
171 : : //choose the scheme for edge mesher
172 [ + - + - : 287 : switch(scheme)
+ - - ]
173 : : {
174 : : case EQUAL://equal meshing for edges
175 [ + - ]: 264 : EqualMeshing(me, num_edges, coords);
176 : 264 : break;
177 : : case BIAS://bias meshing for edges
178 [ # # ]: 0 : BiasMeshing(me, num_edges, coords);
179 : 0 : break;
180 : : case DUAL://dual bias meshing for edges
181 [ + - ]: 4 : DualBiasMeshing(me, num_edges, coords);
182 : 4 : break;
183 : : case CURVATURE://curvature-based meshing for edges
184 [ # # ]: 0 : CurvatureMeshing(me, num_edges, coords);
185 : 0 : break;
186 : : case VARIABLE: // use a var size from sizing function
187 [ + - ]: 19 : VariableMeshing(me, num_edges, coords);
188 : 19 : break;
189 : : case EQUIGNOMONIC: // used to generate HOMME type meshes on a sphere
190 [ # # ]: 0 : EquiAngleGnomonic(me, num_edges, coords);
191 : 0 : break;
192 : : default:
193 : 0 : break;
194 : : }
195 : : //the variable nodes should be resized, node size may be changed in the different scheme for edge meshing
196 [ + - ]: 287 : me->mesh_intervals(num_edges);
197 [ + - ]: 287 : nodes.resize(num_edges+1);
198 : :
199 : : //move the other nodes to the end of nodes' list
200 [ + + ][ + - ]: 287 : if (periodic) nodes[num_edges] = nodes[0];
[ + - ]
201 [ + - ][ + - ]: 271 : else nodes[num_edges] = nodes[1];
202 : :
203 : : //create the vertices' entities on the edge
204 [ + + ]: 287 : if (num_edges > 1) {
205 [ + - ][ + - ]: 239 : rval = mk_core()->moab_instance()->create_vertices(&coords[3], num_edges - 1, mit -> second);
[ + - ][ + - ]
[ + - ]
206 [ - + ][ # # ]: 239 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
207 : : }
208 : :
209 : : //distribute the nodes into vector
210 : 287 : int j = 1;
211 [ + - ][ + - ]: 2489 : for (moab::Range::iterator rit = mit -> second.begin(); rit != mit -> second.end(); rit++)
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ]
212 [ + - ][ + - ]: 2202 : nodes[j++] = *rit;
213 : :
214 : : //get the query interface, which we will use to create the edges directly
215 : : moab::ReadUtilIface *iface;
216 [ + - ][ + - ]: 287 : rval = mk_core() -> moab_instance() -> query_interface(iface);
[ + - ]
217 [ - + ][ # # ]: 287 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
218 : :
219 : : //create the edges, get a direct ptr to connectivity
220 : : moab::EntityHandle starth, *connect, *tmp_connect;
221 [ + - ]: 287 : rval = iface -> get_element_connect(num_edges, 2, moab::MBEDGE, 1, starth, connect);
222 [ - + ][ # # ]: 287 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
223 : :
224 : : //add edges to range for the MESelection
225 [ + - ][ + - ]: 287 : mit -> second.insert(starth, starth + num_edges - 1);
226 : :
227 : : //now set the connectivity array from the nodes
228 [ + - ]: 287 : tmp_connect = &nodes[0];
229 [ + + ]: 2776 : for (int i = 0; i < num_edges; i++)
230 : : {
231 : 2489 : connect[0] = tmp_connect[0];
232 : 2489 : connect[1] = tmp_connect[1];
233 : :
234 : : //increment connectivity ptr by 2 to go to next edge
235 : 2489 : connect += 2;
236 : :
237 : : //increment tmp_connect by 1, to go to next node
238 : 2489 : tmp_connect++;
239 : : }
240 : :
241 : : // ok, we are done, commit to ME
242 [ + - ][ + - ]: 287 : me->commit_mesh(mit->second, COMPLETE_MESH);
243 : :
244 : :
245 : 61 : }
246 : :
247 : :
248 : 61 : }
249 : :
250 : : //---------------------------------------------------------------------------//
251 : : // Deconstruction function
252 : 180 : EdgeMesher::~EdgeMesher()
253 : : {
254 : :
255 [ - + ]: 120 : }
256 : :
257 : : //---------------------------------------------------------------------------//
258 : : // Create the mesh for edges with the equal distances
259 : 264 : void EdgeMesher::EqualMeshing(ModelEnt *ent, int num_edges, std::vector<double> &coords)
260 : : {
261 : : double umin, umax, measure;
262 : : (void) measure;
263 : : //get the u range for the edge
264 [ + - ][ + - ]: 264 : iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax);
[ + - ]
265 [ + - ]: 264 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
266 : :
267 [ - + ][ # # ]: 264 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
268 : :
269 : : //get the arc length
270 [ + - ]: 264 : measure = ent -> measure();
271 : :
272 : : double u, du;
273 [ - + ][ # # ]: 264 : if (!num_edges) throw Error(MK_BAD_INPUT, "Trying to mesh edge with zero edges.");
274 : 264 : du = (umax - umin)/(double)num_edges;
275 : :
276 : 264 : u = umin;
277 : : //distribute the nodes with equal distances
278 [ + + ]: 2271 : for (int i = 1; i < num_edges; i++)
279 : : {
280 : 2007 : u = umin + i*du;
281 [ + - ][ + - ]: 2007 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]);
[ + - ][ + - ]
[ + - ][ + - ]
282 [ + - ]: 2007 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
283 : : }
284 : :
285 : 264 : }
286 : :
287 : : //---------------------------------------------------------------------------//
288 : : // Create the mesh for edges based on the curvatures
289 : 0 : void EdgeMesher::CurvatureMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords)
290 : : {
291 : : double umin, umax, measure;
292 : : (void) measure;
293 : : //store the initial edge size, the edge size may be changed during meshing
294 : 0 : int initial_num_edges = num_edges;
295 : :
296 : : //get the u range for the edge
297 [ # # ][ # # ]: 0 : iGeom::Error gerr = ent->igeom_instance() ->getEntURange(ent->geom_handle(), umin, umax);
[ # # ]
298 [ # # ]: 0 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
299 : :
300 [ # # ][ # # ]: 0 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
301 : :
302 : : //get the arc length
303 [ # # ]: 0 : measure = ent -> measure();
304 : :
305 : 0 : int index = 0;
306 : : double u, du, uMid;
307 : 0 : du = (umax - umin)/(double)num_edges;
308 : :
309 [ # # ]: 0 : std::vector<double> NodeCoordinates;
310 [ # # ]: 0 : std::vector<double> TempNode;
311 [ # # ]: 0 : std::vector<double> URecord; //record the value of U
312 : :
313 : : Point3D pts0, pts1, ptsMid;
314 : : double tmp[3];
315 : :
316 [ # # ]: 0 : NodeCoordinates.resize(3*(num_edges + 1));
317 : :
318 [ # # ]: 0 : TempNode.resize(3*1);
319 [ # # ]: 0 : URecord.resize(1);
320 : :
321 [ # # ][ # # ]: 0 : gerr = ent -> igeom_instance() -> getEntUtoXYZ(ent->geom_handle(), umin, TempNode[0], TempNode[1], TempNode[2]);
[ # # ][ # # ]
[ # # ][ # # ]
322 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge");
323 : :
324 [ # # ][ # # ]: 0 : NodeCoordinates[3*0] = TempNode[0];
325 [ # # ][ # # ]: 0 : NodeCoordinates[3*0+1] = TempNode[1];
326 [ # # ][ # # ]: 0 : NodeCoordinates[3*0+2] = TempNode[2];
327 : :
328 [ # # ]: 0 : URecord[0] = umin;
329 : :
330 : 0 : u = umin;
331 : : //distribute the mesh nodes on the edge based on the curvature
332 [ # # ]: 0 : for (int i = 1; i < num_edges; i++)
333 : : {
334 : : //first distribute the nodes evenly
335 : 0 : u = umin + i*du;
336 [ # # ][ # # ]: 0 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, NodeCoordinates[3*i], NodeCoordinates[3*i+1], NodeCoordinates[3*i+2]);
[ # # ][ # # ]
[ # # ][ # # ]
337 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
338 : :
339 : : //store the two adjacent mesh nodes on the edge
340 [ # # ]: 0 : pts0.px = NodeCoordinates[3*(i-1)];
341 [ # # ]: 0 : pts0.py = NodeCoordinates[3*(i-1)+1];
342 [ # # ]: 0 : pts0.pz = NodeCoordinates[3*(i-1)+2];
343 : :
344 [ # # ]: 0 : pts1.px = NodeCoordinates[3*i];
345 [ # # ]: 0 : pts1.py = NodeCoordinates[3*i+1];
346 [ # # ]: 0 : pts1.pz = NodeCoordinates[3*i+2];
347 : :
348 : : //get the coordinates for mid point between two adjacent mesh nodes on the edge
349 : 0 : uMid = (u-du+u)/2;
350 [ # # ][ # # ]: 0 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uMid, tmp[0], tmp[1], tmp[2]);
[ # # ]
351 : 0 : ptsMid.px = tmp[0];
352 : 0 : ptsMid.py = tmp[1];
353 : 0 : ptsMid.pz = tmp[2];
354 : :
355 : : //calculate the error and check whether it requires further refinement based on the curvature
356 [ # # ][ # # ]: 0 : if (!ErrorCalculate(ent, pts0, pts1, ptsMid))
357 : : {
358 [ # # ]: 0 : DivideIntoMore(ent, pts0, ptsMid, pts1, u-du, u, uMid, index, TempNode, URecord);
359 : : }
360 : :
361 : : // add the other end node to the array, get the next two adjacent mesh nodes
362 : 0 : index++;
363 [ # # ]: 0 : TempNode.resize(3*(index + 1));
364 [ # # ]: 0 : URecord.resize(index + 1);
365 [ # # ]: 0 : TempNode[3*index] = pts1.px;
366 [ # # ]: 0 : TempNode[3*index + 1] = pts1.py;
367 [ # # ]: 0 : TempNode[3*index + 2] = pts1.pz;
368 : :
369 [ # # ]: 0 : URecord[index] = u;
370 : : }
371 : :
372 : : //sorting the parametrical coordinate data based on the value of u
373 [ # # ]: 0 : assert(TempNode.size()== (3*URecord.size()) );
374 : :
375 [ # # ]: 0 : QuickSorting(TempNode, URecord, URecord.size());
376 : 0 : num_edges = URecord.size() - 1;
377 : :
378 : : //resize the variable coords
379 [ # # ]: 0 : coords.resize(3*(num_edges+1));
380 : :
381 : : //move the other end node to the endmost of the list
382 [ # # ]: 0 : for (int i = 0; i < 3; i++)
383 [ # # ][ # # ]: 0 : coords[3*num_edges+i] = coords[3*initial_num_edges+i];
384 : :
385 : : //return the mesh nodes
386 [ # # ]: 0 : for (int i = 1; i < num_edges; i++)
387 : : {
388 [ # # ][ # # ]: 0 : coords[3*i] = TempNode[3*i];
389 [ # # ][ # # ]: 0 : coords[3*i+1] = TempNode[3*i+1];
390 [ # # ][ # # ]: 0 : coords[3*i+2] = TempNode[3*i+2];
391 : 0 : }
392 : :
393 : 0 : }
394 : :
395 : : //---------------------------------------------------------------------------//
396 : : // Create the mesh for edges with dual bias distances
397 : 4 : void EdgeMesher::DualBiasMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords)
398 : : {
399 : : double umin, umax, measure;
400 : :
401 : : //get the u range for the edge
402 [ + - ][ + - ]: 4 : iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax);
[ + - ]
403 [ + - ]: 4 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
404 : :
405 [ - + ][ # # ]: 4 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
406 : :
407 : : //get the arc length
408 [ + - ]: 4 : measure = ent -> measure();
409 : :
410 : : double u, L0, dist, u0, u1;
411 : :
412 : : //if the node # is odd, node # will increase by 1
413 [ - + ]: 4 : if ((num_edges%2)!=0)
414 : : {
415 : 0 : num_edges++;
416 [ # # ]: 0 : coords.resize(3*(num_edges+1));
417 : : //move the other end node's position because the variable coords has been resized.
418 [ # # ]: 0 : for (int k = 0; k < 3; k++)
419 [ # # ][ # # ]: 0 : coords[3*num_edges + k] = coords[3*num_edges + k - 3];
420 : : }
421 : :
422 : : //default bias ratio is 1.2
423 : 4 : double q = ratio;//
424 : :
425 : : //get the distance between the first two nodes
426 [ + - ]: 4 : L0 = 0.5 * measure * (1-q) / (1 - pow(q, num_edges/2));
427 : :
428 : :
429 : 4 : u0 = umin;
430 : 4 : u1 = umax;
431 : : //distribute the mesh nodes on the edge with dual-bias distances
432 [ + + ]: 24 : for (int i = 1; i < (num_edges/2 + 1); i++)
433 : : {
434 : : //distribute the mesh nodes on the one side
435 [ + - ]: 20 : dist = L0*pow(q, i-1);
436 : 20 : u = u0 + (umax - umin) * dist/measure;
437 [ + - ]: 20 : u = getUCoord(ent, u0, dist, u, umin, umax);
438 : 20 : u0 = u;
439 [ + - ][ + - ]: 20 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]);
[ + - ][ + - ]
[ + - ][ + - ]
440 [ + - ]: 20 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
441 : :
442 : : //distribute the mesh nodes on the other side
443 [ + + ]: 20 : if (i < num_edges/2)
444 : : {
445 : 16 : u = u1 - (umax-umin) * dist / measure;
446 [ + - ]: 16 : u = getUCoord(ent, u1, dist, u, umin, umax);
447 : 16 : u1 = u;
448 [ + - ][ + - ]: 16 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*(num_edges-i)], coords[3*(num_edges-i)+1], coords[3*(num_edges-i)+2]);
[ + - ][ + - ]
[ + - ][ + - ]
449 [ + - ]: 16 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge");
450 : : }
451 : :
452 : : }
453 : :
454 : 4 : }
455 : :
456 : : //---------------------------------------------------------------------------//
457 : : // Create the mesh for edges with bias distances
458 : 0 : void EdgeMesher::BiasMeshing(ModelEnt *ent, int num_edges, std::vector<double> &coords)
459 : : {
460 : : double umin, umax, measure;
461 : :
462 : : //get the u range for the edge
463 [ # # ][ # # ]: 0 : iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax);
[ # # ]
464 [ # # ]: 0 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
465 : :
466 [ # # ][ # # ]: 0 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
467 : :
468 : : //get the arc length
469 [ # # ]: 0 : measure = ent -> measure();
470 : :
471 : 0 : double u, L0, dist = 0, u0;
472 : :
473 : : // the default bias ratio 1.2
474 : 0 : double q = ratio;
475 [ # # ]: 0 : L0 = measure * (1-q) / (1 - pow(q, num_edges));
476 : :
477 : :
478 : 0 : u0 = umin;
479 : : //distribute the mesh nodes on the edge with bias distances
480 [ # # ]: 0 : for (int i = 1; i < num_edges; i++)
481 : : {
482 [ # # ]: 0 : dist = L0*pow(q, i-1);
483 : 0 : u = u0 + (umax - umin)*dist/measure;
484 [ # # ]: 0 : u = getUCoord(ent, u0, dist, u, umin, umax);
485 : 0 : u0 = u;
486 [ # # ][ # # ]: 0 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]);
[ # # ][ # # ]
[ # # ][ # # ]
487 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
488 : : }
489 : 0 : }
490 : :
491 : : //---------------------------------------------------------------------------//
492 : : // Mesh Refinement function based on the curvature: if the error is too big, refine the mesh on the edge
493 : 0 : void EdgeMesher::DivideIntoMore(ModelEnt *ent, Point3D p0, Point3D pMid, Point3D p1, double u0, double u1, double uMid, int &index, vector<double> &nodes, vector<double> &URecord)
494 : : {
495 : : //this is a recursive process, the process continues until the error is smaller than what is required.
496 : : //first get two adjacent mesh nodes on the edge and coordinates for mid point between two adjacent nodes.
497 : : //then check the left side and right side whether the error is too big nor not
498 : : double uu0, uu1, uumid, tmp[3];
499 : : Point3D pts0, pts1, ptsMid;
500 : :
501 : 0 : index++;
502 [ # # ]: 0 : nodes.resize(3*(index+1));
503 [ # # ]: 0 : URecord.resize(index+1);
504 [ # # ]: 0 : nodes[3*index] = pMid.px;
505 [ # # ]: 0 : nodes[3*index+1] = pMid.py;
506 [ # # ]: 0 : nodes[3*index+2] = pMid.pz;
507 [ # # ]: 0 : URecord[index] = uMid;
508 : :
509 : : //check the left side
510 : 0 : uu0=u0;
511 : 0 : uu1=uMid;
512 : 0 : uumid=(uu0+uu1)/2;
513 : 0 : pts0=p0;
514 : 0 : pts1=pMid;
515 : :
516 : :
517 [ # # ][ # # ]: 0 : iGeom::Error gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uumid, tmp[0], tmp[1], tmp[2]);
[ # # ]
518 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
519 : 0 : ptsMid.px = tmp[0];
520 : 0 : ptsMid.py = tmp[1];
521 : 0 : ptsMid.pz = tmp[2];
522 : :
523 : : //check the error
524 [ # # ][ # # ]: 0 : if(!ErrorCalculate(ent, pts0, pts1, ptsMid))
525 : : {
526 : : //further refinement
527 [ # # ]: 0 : DivideIntoMore(ent, pts0, ptsMid, pts1, uu0, uu1, uumid, index, nodes, URecord);
528 : : }
529 : :
530 : : //check the right side
531 : 0 : uu0 = uMid;
532 : 0 : uu1=u1;
533 : 0 : uumid=(uu0+uu1)/2;
534 : 0 : pts0=pMid;
535 : 0 : pts1=p1;
536 : : //get the coorindinates for mid point
537 [ # # ][ # # ]: 0 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), uumid, tmp[0], tmp[1], tmp[2]);
[ # # ]
538 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
539 : 0 : ptsMid.px = tmp[0];
540 : 0 : ptsMid.py = tmp[1];
541 : 0 : ptsMid.pz = tmp[2];
542 : :
543 : : //check the error
544 [ # # ][ # # ]: 0 : if(!ErrorCalculate(ent, pts0, pts1, ptsMid))
545 : : {
546 : : //further refinement
547 [ # # ]: 0 : DivideIntoMore(ent, pts0, ptsMid, pts1, uu0, uu1, uumid, index, nodes, URecord);
548 : : }
549 : :
550 : 0 : }
551 : :
552 : : //create the mesh for edges based on variable size from SizingFunction (var)
553 : 19 : void EdgeMesher::VariableMeshing(ModelEnt *ent, int &num_edges, std::vector<double> &coords)
554 : : {
555 : : double umin, umax, measure;
556 : : (void) measure;
557 : : // because of that, keep track of the first node position and last node position
558 : : // first node position does not change, but the last node position do change
559 : : // coords will contain all nodes, including umax in Urecord!
560 : :
561 [ + - ][ + - ]: 19 : SizingFunction *sf = mk_core()->sizing_function(ent->sizing_function_index());
[ + - ]
562 : : //get the u range for the edge
563 [ + - ]: 19 : iGeom::EntityHandle edge = ent->geom_handle();
564 [ + - ][ + - ]: 19 : iGeom::Error gerr = ent->igeom_instance() ->getEntURange(edge, umin, umax);
565 [ + - ]: 19 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
566 : :
567 [ - + ][ # # ]: 19 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
568 : :
569 : : //get the arc length
570 [ + - ]: 19 : measure = ent -> measure();
571 : :
572 : : // start advancing for each edge mesh, from the first point position
573 : 19 : double currentPar = umin;
574 : : double currentPosition[3];
575 [ + - ]: 19 : gerr = ent->igeom_instance() ->getEntUtoXYZ(edge, umin, currentPosition[0],
576 [ + - ]: 19 : currentPosition[1], currentPosition[2] );
577 : :
578 : : double endPoint[3];
579 [ + - ]: 19 : gerr = ent->igeom_instance() ->getEntUtoXYZ(edge, umax, endPoint[0],
580 [ + - ]: 19 : endPoint[1], endPoint[2] );
581 [ + - ]: 19 : Vector<3> endpt(endPoint);
582 : :
583 [ + - ]: 19 : double targetSize = sf->size(currentPosition);
584 : 19 : double startSize = targetSize;
585 : :
586 [ + - ]: 19 : double endSize = sf->size(endPoint);
587 : : // advance u such that the next point is at "currentSize" distance
588 : : // or close to it
589 : : // try first with a u that is coming from the (umax-umin)/number of edges
590 : 19 : double deltaU = (umax-umin)/num_edges;
591 : : //coords.clear(); we do not want to clear, as the first node is still fine
592 [ + - ]: 19 : std::vector<double> URecord; //record the values for u; we may have to adjust all
593 : : // of them accordingly, and even add one more if we have some evenify problems.
594 : : // keep in mind that size is just a suggestion, number of intervals is more important for
595 : : // paver mesher
596 [ + - ]: 19 : Vector<3> pt(currentPosition);
597 : :
598 : : //bool notDone = true;
599 : 19 : double prevU = umin;
600 [ + + ]: 174 : while (currentPar + 1.1*deltaU < umax)
601 : : {
602 : : // do some binary search; actually, better, do Newton-Raphson, which should converge
603 : : // faster
604 : : //
605 : 155 : prevU = currentPar;
606 : 155 : currentPar += deltaU;
607 : : // adjust current par, such that
608 : : double point[3];
609 [ + - ][ + - ]: 155 : gerr=ent->igeom_instance()->getEntUtoXYZ(edge, currentPar, point[0], point[1], point[2] );
610 [ + - ]: 155 : IBERRCHK(gerr, "Trouble getting position at parameter u.");
611 [ + - ]: 155 : Vector<3> ptCandidate(point);
612 [ + - ][ + - ]: 155 : double compSize = length(ptCandidate-pt);
613 : 155 : int nit = 0;
614 : :
615 [ + + ][ + - ]: 303 : while ( (fabs(1.-compSize/targetSize)> 0.02 ) && (nit < 5))// 2% of the target size
616 : : {
617 : : // do Newton iteration
618 : : double tangent[3];
619 [ + - ][ + - ]: 148 : gerr=ent->igeom_instance() ->getEntTgntU(edge, currentPar, tangent[0], tangent[1], tangent[2] );
620 [ + - ]: 148 : IBERRCHK(gerr, "Trouble getting tangent at parameter u.");
621 [ + - ]: 148 : Vector<3> tang(tangent);
622 [ + - ][ + - ]: 148 : double dldu = 1./compSize * ((ptCandidate-pt )%tang);
623 : 148 : nit++;// increase iteration count
624 [ + - ]: 148 : if (dldu!=0.)
625 : : {
626 : 148 : double deu= (targetSize-compSize)/dldu;
627 : 148 : currentPar+=deu;
628 [ - + ]: 148 : if (prevU>currentPar)
629 : : {
630 : 148 : break; // something is wrong
631 : : }
632 [ - + ]: 148 : if (umax < currentPar)
633 : : {
634 : 0 : currentPar = umax;
635 : 0 : break;
636 : : }
637 [ + - ][ + - ]: 148 : ent->igeom_instance()->getEntUtoXYZ(edge, currentPar, point[0], point[1], point[2]);
638 [ + - ]: 148 : Vector<3> newPt(point);
639 [ + - ][ + - ]: 148 : compSize = length(newPt-pt);
640 : 148 : ptCandidate = newPt;
641 : : }
642 : :
643 : : }
644 : : // we have found an acceptable point/param
645 [ + - ]: 155 : URecord.push_back(currentPar);
646 : 155 : deltaU = currentPar-prevU;// should be greater than 0
647 : 155 : pt = ptCandidate;
648 [ + - ][ + - ]: 155 : targetSize = sf->size(pt.data());// a new target size, at the current point
649 : :
650 : :
651 : :
652 : : }
653 : : // when we are out of here, we need to adjust the URecords, to be more nicely
654 : : // distributed; also, look at the evenify again
655 : 19 : int sizeU = (int)URecord.size();
656 [ + + ][ + - ]: 19 : if ((sizeU%2==0) && ent->constrain_even() )
[ + + ][ + + ]
657 : : {
658 : : // add one more
659 [ - + ]: 4 : if (sizeU==0)
660 : : {
661 : : // just add one in the middle, and call it done
662 [ # # ]: 0 : URecord.push_back( (umin+umax)/2);
663 : : }
664 : : else
665 : : {
666 : : //at least 2 (maybe 4)
667 [ + - ][ + - ]: 4 : double lastDelta = URecord[sizeU-1]-URecord[sizeU-2];
668 [ + - ][ + - ]: 4 : URecord.push_back(URecord[sizeU-1]+lastDelta );
669 : : }
670 : : }
671 : : // now, we have to redistribute, such as the last 2 deltas are about the same
672 : : // so, we should have after a little work,
673 : : // umin, umin+c*(URecord[0]-umin), ... umin+c*(URecord[size-1]-umin), umax
674 : : // what we have now is
675 : : // umin, URecord[0], ... ,URecord[size-1], and umax could be even outside or inside
676 : : // keep the last sizes equal
677 : : // umin+c(UR[size-2]-umin) + umax = 2*( umin+c*(UR[size-1]-umin))
678 : : // c ( 2*(UR[size-1]-umin) -(UR[size-2]-umin) ) = umax - umin
679 : : // c ( 2*UR[size-1] - UR[size-2] - umin ) = umax - umin
680 : : // c = (umax-umin)/( 2*UR[size-1] - UR[size-2] - umin)
681 : 19 : sizeU = (int)URecord.size();// it may be bigger by one than the last time
682 [ + - ]: 19 : if (sizeU == 0)
683 : : {
684 : : // nothing to do, only one edge to generate
685 : : }
686 [ - + ]: 19 : else if (sizeU == 1)
687 : : {
688 : : // put it according to the sizes at ends, and assume a nice variation for u
689 : : // (u-umin) / (umax-u) = startSize / endSize
690 : : // (u-umin)*endSize = (umax-u) * startSize
691 : : // u(endSize+startSize)=(umax*startSize+umin*endSize)
692 [ # # ]: 0 : URecord[0] = (umax*startSize+umin*endSize)/(startSize+endSize);
693 : :
694 : : }
695 : : else // sizeU>=2, so we can spread the param a little more, assuming nice
696 : : // uniform mapping
697 : : {
698 [ + - ][ + - ]: 19 : double c = (umax-umin)/( 2*URecord[sizeU-1] - URecord[sizeU-2] - umin);
699 [ + + ]: 178 : for (int i=0; i<sizeU; i++)
700 [ + - ][ + - ]: 159 : URecord[i] = umin + c*(URecord[i] -umin);// some spreading out
701 : : }
702 : : // now, we can finally get the points for each U, U's should be spread out nicely
703 [ + - ]: 19 : URecord.push_back(umax); // just add the last u, for the end point
704 : : //
705 : 19 : sizeU = (int) URecord.size(); // new size, after pushing the last u, umax
706 : 19 : num_edges = sizeU;// this is the new number of edges; the last one will be the end point
707 : : // of the edge, corresponding to umax
708 [ + - ]: 19 : coords.resize(3*sizeU+3);
709 : : // we already know that at i=0 is the first node, start vertex of edge
710 : : // the rest will be computed from u
711 : : // even the last one, which is an overkill
712 [ + + ]: 197 : for (int i = 1; i <= num_edges; i++)
713 : : {
714 [ + - ]: 178 : double u = URecord[i-1];
715 [ + - ][ + - ]: 178 : gerr = ent->igeom_instance()->getEntUtoXYZ(edge, u, coords[3*i], coords[3*i+1], coords[3*i+2]);
[ + - ][ + - ]
[ + - ]
716 [ + - ]: 178 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
717 : : }
718 : 19 : return;
719 : :
720 : : }
721 : : // number of edges is input here
722 : : // equal angles are formed at the center of the sphere/cube mesh
723 : : // it is close to bias meshing, but not quite
724 : 0 : void EdgeMesher::EquiAngleGnomonic(ModelEnt *ent, int num_edges, std::vector<double> &coords)
725 : : {
726 : 0 : const double MY_PI=3.14159265;
727 : 0 : double deltaAngle=MY_PI/num_edges/2;
728 : : // double length=ent->measure();// this is an edge
729 : : double umin, umax;
730 : :
731 : : //get the u range for the edge
732 [ # # ][ # # ]: 0 : iGeom::Error gerr = ent->igeom_instance()->getEntURange(ent->geom_handle(), umin, umax);
[ # # ]
733 [ # # ]: 0 : IBERRCHK(gerr, "Trouble get parameter range for edge.");
734 : :
735 [ # # ][ # # ]: 0 : if (umin == umax) throw Error(MK_BAD_GEOMETRIC_EVALUATION, "Edge evaluated to some parameter umax and umin.");
736 : :
737 : : // consider that the parametrization is very linear
738 : : // most of the time u will be from 0 to length of edge, for a cube
739 : 0 : double deltau = umax - umin;
740 : :
741 : 0 : double u = umin;// u will get different values,
742 : : // start at u
743 [ # # ]: 0 : for (int i = 1; i < num_edges; i++)
744 : : {
745 : 0 : double betak=i*deltaAngle;
746 : 0 : double alfak = MY_PI/4-betak;
747 : 0 : double tang_alfak = tan(alfak);
748 : 0 : u = umin+deltau/2*(1-tang_alfak);
749 : :
750 [ # # ][ # # ]: 0 : gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, coords[3*i], coords[3*i+1], coords[3*i+2]);
[ # # ][ # # ]
[ # # ][ # # ]
751 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
752 : : }
753 : 0 : return;
754 : : }
755 : : //---------------------------------------------------------------------------//
756 : : // Rapid sorting the mesh nodes on the edge based on the parametric coordinates. This is a recursive
757 : : // process
758 : 0 : void EdgeMesher::RapidSorting(vector<double> &nodes, vector<double> &URecord, int left, int right)
759 : : {
760 : : int i, j;
761 : : double middle, iTemp;
762 : : Point3D TempData;
763 : :
764 [ # # ]: 0 : middle=URecord[(left+right)/2];
765 : 0 : i=left;
766 : 0 : j=right;
767 : :
768 [ # # ]: 0 : do
769 : : {
770 : : //search the values which are greater than the middle value from the left side
771 [ # # ][ # # ]: 0 : while((URecord[i] < middle)&&(i<right))
[ # # ][ # # ]
772 : : {
773 : 0 : i++;
774 : : }
775 : : //search the values which are greater than the middle value from the right side
776 [ # # ][ # # ]: 0 : while((URecord[j] > middle)&&(j > left))
[ # # ][ # # ]
777 : : {
778 : 0 : j--;
779 : : }
780 [ # # ]: 0 : if (i<=j)//find a pair of values
781 : : {
782 [ # # ]: 0 : iTemp = URecord[i];
783 [ # # ][ # # ]: 0 : URecord[i] = URecord[j];
784 [ # # ]: 0 : URecord[j]=iTemp;
785 : :
786 : :
787 [ # # ]: 0 : TempData.px = nodes[3*i];
788 [ # # ]: 0 : TempData.py = nodes[3*i+1];
789 [ # # ]: 0 : TempData.pz = nodes[3*i+2];
790 : :
791 [ # # ][ # # ]: 0 : nodes[3*i] = nodes[3*j];
792 [ # # ][ # # ]: 0 : nodes[3*i+1] = nodes[3*j+1];
793 [ # # ][ # # ]: 0 : nodes[3*i+2] = nodes[3*j+2];
794 [ # # ]: 0 : nodes[3*j] = TempData.px;
795 [ # # ]: 0 : nodes[3*j+1] = TempData.py;
796 [ # # ]: 0 : nodes[3*j+2] = TempData.pz;
797 : :
798 : 0 : i++;
799 : 0 : j--;
800 : : }
801 : : }while(i<=j);
802 [ # # ]: 0 : if (left < j)
803 [ # # ]: 0 : RapidSorting(nodes, URecord, left, j);
804 [ # # ]: 0 : if (right > i)
805 [ # # ]: 0 : RapidSorting(nodes, URecord, i, right);
806 : 0 : }
807 : :
808 : : //---------------------------------------------------------------------------//
809 : : // Quick Sorting: this function comes together with the function RapidSorting
810 : 0 : void EdgeMesher::QuickSorting(vector<double> &nodes, vector<double> &URecord, int count)
811 : : {
812 : 0 : RapidSorting(nodes, URecord, 0, count-1);
813 : 0 : }
814 : :
815 : :
816 : : //---------------------------------------------------------------------------//
817 : : // return the x, y, z coordinates from the parametric coordinates
818 : 72 : EdgeMesher::Point3D EdgeMesher::getXYZCoords(ModelEnt *ent, double u) const
819 : : {
820 : : Point3D pts3D;
821 : : double xyz[3];
822 : :
823 : :
824 : : //get the coordinates in the physical space
825 [ + - ][ + - ]: 72 : iGeom::Error gerr = ent->igeom_instance()->getEntUtoXYZ(ent->geom_handle(), u, xyz[0], xyz[1], xyz[2]);
[ + - ]
826 [ + - ]: 72 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
827 : :
828 : 72 : pts3D.px = xyz[0];
829 : 72 : pts3D.py = xyz[1];
830 : 72 : pts3D.pz = xyz[2];
831 : 72 : return pts3D;
832 : : }
833 : :
834 : : //---------------------------------------------------------------------------//
835 : : // get the parametric coordinate based on first parametric coordinate ustart and distance in the physical space
836 : 36 : double EdgeMesher::getUCoord(ModelEnt *ent, double ustart, double dist, double uguess, double umin, double umax) const
837 : : {
838 : :
839 [ + - ]: 36 : Point3D p0 = getXYZCoords(ent, ustart);
840 [ + - ]: 36 : Point3D p1 = getXYZCoords(ent, uguess);
841 : :
842 : :
843 : 36 : double dx, dy, dz, dl, u=uguess;
844 : 36 : double tol = 1.0E-7;
845 : 36 : int test=0;
846 : :
847 : 36 : int ntrials=0;
848 : : while(1)
849 : : {
850 : 36 : dx = p1.px - p0.px;
851 : 36 : dy = p1.py - p0.py;
852 : 36 : dz = p1.pz - p0.pz;
853 : 36 : dl = sqrt(dx * dx + dy * dy + dz * dz);
854 [ + - ]: 36 : if ( fabs(dl-dist) < tol) break;
855 : :
856 : 0 : u = ustart + (u - ustart) * (dist/dl);
857 [ # # ]: 0 : if (u > umax)
858 : : {
859 : 0 : u=umax;
860 : 0 : test++;
861 [ # # ]: 0 : if (test>10) break;
862 : : }
863 [ # # ]: 0 : if (u < umin)
864 : : {
865 : 0 : u=umin;
866 : 0 : test++;
867 [ # # ]: 0 : if (test>10) break;
868 : : }
869 [ # # ]: 0 : p1 = getXYZCoords(ent, u);
870 : :
871 : :
872 [ # # ]: 0 : if (ntrials++ == 100000)
873 : : {
874 [ # # ][ # # ]: 0 : cout << " Warning: Searching for U failed " << endl;
875 : : }
876 : : }
877 : 36 : uguess = u;
878 : 36 : return uguess;
879 : : }
880 : :
881 : : //---------------------------------------------------------------------------//
882 : : // calculate the error: the distance between the mid point of two adjacent
883 : : // mesh nodes (on the mesh line segments) and mid point on the edge
884 : 0 : bool EdgeMesher::ErrorCalculate(ModelEnt *ent, Point3D p0, Point3D p1, Point3D pMid)
885 : : {
886 : : double lengtha, lengthb, lengthc;
887 : : double deltax, deltay, deltaz;
888 : 0 : double angle, error, tol=1.0E-3, H;
889 : : double cvtr_ijk[3], curvature;
890 : : bool result;
891 : :
892 : : //calculate the distance between the first mesh node and mid point on the edge
893 : 0 : deltax = pMid.px-p0.px;
894 : 0 : deltay = pMid.py-p0.py;
895 : 0 : deltaz = pMid.pz-p0.pz;
896 : 0 : lengtha = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz);
897 : :
898 : : //calculate the distance between two adjacent mesh nodes
899 : 0 : deltax = p1.px-p0.px;
900 : 0 : deltay = p1.py-p0.py;
901 : 0 : deltaz = p1.pz-p0.pz;
902 : 0 : lengthb = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz);
903 : :
904 : : //calculate the distance between the second mesh node and mid point on the edge
905 : 0 : deltax = pMid.px-p1.px;
906 : 0 : deltay = pMid.py-p1.py;
907 : 0 : deltaz = pMid.pz-p1.pz;
908 : 0 : lengthc = sqrt(deltax*deltax + deltay*deltay + deltaz*deltaz);
909 : :
910 : : //calculate the angle
911 : 0 : angle = acos((lengtha*lengtha + lengthb*lengthb - lengthc*lengthc)/(2*lengtha*lengthb));
912 : 0 : H = fabs(lengtha*sin(angle));
913 : :
914 : : //calculate the curvature
915 [ # # ][ # # ]: 0 : iGeom::Error gerr = ent->igeom_instance()->getEgCvtrXYZ(ent->geom_handle(), pMid.px, pMid.py, pMid.pz, cvtr_ijk[0], cvtr_ijk[1], cvtr_ijk[2]);
[ # # ]
916 [ # # ]: 0 : IBERRCHK(gerr, "Trouble getting U from XYZ along the edge.");
917 : 0 : curvature = sqrt(cvtr_ijk[0]*cvtr_ijk[0]+cvtr_ijk[1]*cvtr_ijk[1]+cvtr_ijk[2]*cvtr_ijk[2]);
918 : 0 : error= H*curvature;
919 : :
920 : :
921 [ # # ]: 0 : if (error > tol)
922 : 0 : result = false;
923 : : else
924 : 0 : result = true;
925 : 0 : return result;
926 : : }
927 : :
928 [ + - ][ + - ]: 156 : }
929 : :
|
