1 #include <iostream>
2 #include <fstream>
3 #include <algorithm>
4 #include <iomanip>
5 #include <cassert>
6 #include <limits>
7 #include "moab/OrientedBoxTreeTool.hpp"
8 #include "SmoothFace.hpp"
9
10 #define GEOMETRY_RESABS 1.e-6
11 #define mbsqr(a) ((a)*(a))
12 #define mbcube(a) (mbsqr(a) * (a))
13 #define mbquart(a) (mbsqr(a) * mbsqr(a))
14
15 namespace moab {
16
within_tolerance(CartVect & p1,CartVect & p2,const double & tolerance)17 bool within_tolerance(CartVect & p1, CartVect & p2, const double & tolerance)
18 {
19 if ((fabs(p1[0] - p2[0]) < tolerance) && (fabs(p1[1] - p2[1]) < tolerance)
20 && (fabs(p1[2] - p2[2]) < tolerance))
21 return true;
22 return false;
23 }
numAdjTriInSet(Interface * mb,EntityHandle startEdge,EntityHandle set)24 int numAdjTriInSet(Interface * mb, EntityHandle startEdge, EntityHandle set)
25 {
26 std::vector<EntityHandle> adjTri;
27 mb->get_adjacencies(&startEdge, 1, 2, false, adjTri, Interface::UNION);
28 // count how many are in the set
29 int nbInSet = 0;
30 for (size_t i = 0; i < adjTri.size(); i++)
31 {
32 EntityHandle tri = adjTri[i];
33 if (mb->contains_entities(set, &tri, 1))
34 nbInSet++;
35 }
36 return nbInSet;
37 }
38
39 bool debug_surf_eval1 = false;
40
SmoothFace(Interface * mb,EntityHandle surface_set,GeomTopoTool * gTool)41 SmoothFace::SmoothFace(Interface * mb, EntityHandle surface_set,
42 GeomTopoTool * gTool) :
43 _markTag(0), _gradientTag(0), _tangentsTag(0), _edgeCtrlTag(0),
44 _facetCtrlTag(0), _facetEdgeCtrlTag(0), _planeTag(0),
45 _mb(mb), _set(surface_set), _my_geomTopoTool(gTool), _obb_root(0), _evaluationsCounter(0)
46 {
47 //_smooth_face = NULL;
48 //_mbOut->create_meshset(MESHSET_SET, _oSet); //will contain the
49 // get also the obb_root
50 if (_my_geomTopoTool)
51 {
52 _my_geomTopoTool->get_root(this->_set, _obb_root);
53 if (debug_surf_eval1)
54 _my_geomTopoTool->obb_tree()->stats(_obb_root, std::cout);
55 }
56 }
57
~SmoothFace()58 SmoothFace::~SmoothFace()
59 {
60 }
61
area()62 double SmoothFace::area()
63 {
64 // find the area of this entity
65 //assert(_smooth_face);
66 //double area1 = _smooth_face->area();
67 double totArea = 0.;
68 for (Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it)
69 {
70 EntityHandle tria = *it;
71 const EntityHandle * conn3;
72 int nnodes;
73 _mb->get_connectivity(tria, conn3, nnodes);
74 //
75 //double coords[9]; // store the coordinates for the nodes
76 //_mb->get_coords(conn3, 3, coords);
77 CartVect p[3];
78 _mb->get_coords(conn3, 3, (double*) &p[0]);
79 // need to compute the angles
80 // compute angles and the normal
81 //CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
82
83 CartVect AB(p[1] - p[0]);//(p2 - p1);
84 CartVect BC(p[2] - p[1]);//(p3 - p2);
85 CartVect normal = AB * BC;
86 totArea += normal.length() / 2;
87 }
88 return totArea;
89 }
90 // these tags will be collected for deletion
append_smooth_tags(std::vector<Tag> & smoothTags)91 void SmoothFace::append_smooth_tags(std::vector<Tag> & smoothTags)
92 {
93 // these are created locally, for each smooth face
94 smoothTags.push_back(_gradientTag);
95 smoothTags.push_back(_planeTag);
96
97 }
bounding_box(double box_min[3],double box_max[3])98 void SmoothFace::bounding_box(double box_min[3], double box_max[3])
99 {
100
101 for (int i = 0; i < 3; i++)
102 {
103 box_min[i] = _minim[i];
104 box_max[i] = _maxim[i];
105 }
106 // _minim, _maxim
107 }
108
move_to_surface(double & x,double & y,double & z)109 void SmoothFace::move_to_surface(double& x, double& y, double& z)
110 {
111
112 CartVect loc2(x, y, z);
113 bool trim = false;// is it needed?
114 bool outside = true;
115 CartVect closestPoint;
116
117 ErrorCode rval = project_to_facets_main(loc2, trim, outside, &closestPoint,
118 NULL);
119 if (MB_SUCCESS != rval) return;
120 assert(rval==MB_SUCCESS);
121 x = closestPoint[0];
122 y = closestPoint[1];
123 z = closestPoint[2];
124
125 }
126 /*
127 void SmoothFace::move_to_surface(double& x, double& y, double& z,
128 double& u_guess, double& v_guess) {
129 if (debug_surf_eval1) {
130 std::cout << "move_to_surface called." << std::endl;
131 }
132 }*/
133
normal_at(double x,double y,double z,double & nx,double & ny,double & nz)134 bool SmoothFace::normal_at(double x, double y, double z, double& nx,
135 double& ny, double& nz)
136 {
137
138 CartVect loc2(x, y, z);
139
140 bool trim = false;// is it needed?
141 bool outside = true;
142 //CartVect closestPoint;// not needed
143 // not interested in normal
144 CartVect normal;
145 ErrorCode rval = project_to_facets_main(loc2, trim, outside, NULL, &normal);
146 if (MB_SUCCESS != rval) return false;
147 assert(rval==MB_SUCCESS);
148 nx = normal[0];
149 ny = normal[1];
150 nz = normal[2];
151
152 return true;
153 }
154
compute_control_points_on_edges(double min_dot,Tag edgeCtrlTag,Tag markTag)155 ErrorCode SmoothFace::compute_control_points_on_edges(double min_dot,
156 Tag edgeCtrlTag, Tag markTag)
157 {
158
159 _edgeCtrlTag = edgeCtrlTag;
160 _markTag = markTag;
161
162 // now, compute control points for all edges that are not marked already (they are no on the boundary!)
163 for (Range::iterator it = _edges.begin(); it != _edges.end(); ++it)
164 {
165 EntityHandle edg = *it;
166 // is the edge marked? already computed
167 unsigned char tagVal = 0;
168 _mb->tag_get_data(_markTag, &edg, 1, &tagVal);
169 if (tagVal)
170 continue;
171 //double min_dot;
172 init_bezier_edge(edg, min_dot);
173 tagVal = 1;
174 _mb->tag_set_data(_markTag, &edg, 1, &tagVal);
175 }
176 return MB_SUCCESS;
177 }
178
init_gradient()179 int SmoothFace::init_gradient()
180 {
181 // first, create a Tag for gradient (or normal)
182 // loop over all triangles in set, and modify the normals according to the angle as weight
183 // get all the edges from this subset
184 if (NULL == _mb)
185 return 1; //fail
186 _triangles.clear();
187 ErrorCode rval = _mb->get_entities_by_type(_set, MBTRI, _triangles);
188 if (MB_SUCCESS != rval)
189 return 1;
190 // get a new range of edges, and decide the loops from here
191 _edges.clear();
192 rval = _mb-> get_adjacencies(_triangles, 1, true, _edges, Interface::UNION);
193 assert(rval == MB_SUCCESS);
194
195 rval = _mb->get_adjacencies(_triangles, 0, false, _nodes, Interface::UNION);
196 assert(rval == MB_SUCCESS);
197
198 // initialize bounding box limits
199 CartVect vert1;
200 EntityHandle v1 = _nodes[0];// first vertex
201 _mb->get_coords(&v1, 1, (double*) &vert1);
202 _minim = vert1;
203 _maxim = vert1;
204
205 double defNormal[] = { 0., 0., 0. };
206 // look for a tag name here that is definitely unique. We do not want the tags to interfere with each other
207 // this normal will be for each node in a face
208 // some nodes have multiple normals, if they are at the feature edges
209 unsigned long setId = _mb->id_from_handle(_set);
210 char name[50] = { 0 };
211 sprintf(name, "GRADIENT%lu", setId);// name should be something like GRADIENT29, where 29 is the set ID of the face
212 rval = _mb->tag_get_handle(name, 3, MB_TYPE_DOUBLE, _gradientTag,
213 MB_TAG_DENSE|MB_TAG_CREAT, &defNormal);
214 assert(rval == MB_SUCCESS);
215
216 double defPlane[4] = { 0., 0., 1., 0. };
217 // also define a plane tag ; this will be for each triangle
218 char namePlaneTag[50] = { 0 };
219 sprintf(namePlaneTag, "PLANE%lu", setId);
220 rval = _mb->tag_get_handle("PLANE", 4, MB_TYPE_DOUBLE, _planeTag,
221 MB_TAG_DENSE|MB_TAG_CREAT, &defPlane);
222 assert(rval == MB_SUCCESS);
223 // the fourth double is for weight, accumulated at each vertex so far
224 // maybe not needed in the end
225 for (Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it)
226 {
227 EntityHandle tria = *it;
228 const EntityHandle * conn3;
229 int nnodes;
230 _mb->get_connectivity(tria, conn3, nnodes);
231 if (nnodes != 3)
232 return 1; // error
233 //double coords[9]; // store the coordinates for the nodes
234 //_mb->get_coords(conn3, 3, coords);
235 CartVect p[3];
236 _mb->get_coords(conn3, 3, (double*) &p[0]);
237 // need to compute the angles
238 // compute angles and the normal
239 //CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
240
241 CartVect AB(p[1] - p[0]);//(p2 - p1);
242 CartVect BC(p[2] - p[1]);//(p3 - p2);
243 CartVect CA(p[0] - p[2]);//(p1 - p3);
244 double a[3];
245 a[1] = angle(AB, -BC); // angle at B (p2), etc.
246 a[2] = angle(BC, -CA);
247 a[0] = angle(CA, -AB);
248 CartVect normal = -AB * CA;
249 normal.normalize();
250 double plane[4];
251
252 const double * coordNormal = normal.array();
253
254 plane[0] = coordNormal[0];
255 plane[1] = coordNormal[1];
256 plane[2] = coordNormal[2];
257 plane[3] = -normal % p[0]; // dot product
258 //set the plane
259 rval = _mb->tag_set_data(_planeTag, &tria, 1, plane);
260 assert(rval == MB_SUCCESS);
261
262 // add to each vertex the tag value the normal multiplied by the angle
263 double values[9];
264
265 _mb->tag_get_data(_gradientTag, conn3, 3, values);
266 for (int i = 0; i < 3; i++)
267 {
268 //values[4*i]+=a[i]; // first is the weight, which we do not really need
269 values[3 * i + 0] += a[i] * coordNormal[0];
270 values[3 * i + 1] += a[i] * coordNormal[1];
271 values[3 * i + 2] += a[i] * coordNormal[2];
272 }
273
274 // reset those values
275 _mb->tag_set_data(_gradientTag, conn3, 3, values);
276
277 }
278 // normalize the gradients at each node; maybe not needed here?
279 // no, we will do it, it is important
280 int numNodes = _nodes.size();
281 double * normalVal = new double[3 * numNodes];
282 _mb->tag_get_data(_gradientTag, _nodes, normalVal);// get all the normal values at the _nodes
283 for (int i = 0; i < numNodes; i++)
284 {
285 CartVect p1(&normalVal[3 * i]);
286 p1.normalize();
287 p1.get(&normalVal[3 * i]);
288 }
289
290 // reset the normal values after normalization
291 _mb->tag_set_data(_gradientTag, _nodes, normalVal);
292 // print the loops size and some other stuff
293 if (debug_surf_eval1)
294 {
295 std::cout << " normals at " << numNodes << " nodes" << std::endl;
296 int i = 0;
297 for (Range::iterator it = _nodes.begin(); it != _nodes.end(); ++it, i++)
298 {
299 EntityHandle node = *it;
300 std::cout << " Node id " << _mb->id_from_handle(node) << " "
301 << normalVal[3 * i] << " " << normalVal[3 * i + 1] << " "
302 << normalVal[3 * i + 2] << std::endl;
303 }
304 }
305
306 delete[] normalVal;
307
308 return 0;
309 }
310
311 // init bezier edges
312 // start copy
313 //===========================================================================
314 //Function Name: init_bezier_edge
315 //
316 //Member Type: PRIVATE
317 //Description: compute the control points for an edge
318 //===========================================================================
init_bezier_edge(EntityHandle edge,double)319 ErrorCode SmoothFace::init_bezier_edge(EntityHandle edge, double )
320 {
321 // min dot was used for angle here
322 //int stat = 0; // CUBIT_SUCCESS;
323 // all boundaries will be simple, initially
324 // we may complicate them afterwards
325
326 CartVect ctrl_pts[3];
327 int nnodes = 0;
328 const EntityHandle * conn2 = NULL;
329 ErrorCode rval = _mb->get_connectivity(edge, conn2, nnodes);
330 assert(rval == MB_SUCCESS);
331 if (MB_SUCCESS != rval) return rval;
332
333 assert(2 == nnodes);
334 //double coords[6]; // store the coordinates for the nodes
335 CartVect P[2];
336 //ErrorCode rval = _mb->get_coords(conn2, 2, coords);
337 rval = _mb->get_coords(conn2, 2, (double*) &P[0]);
338 assert(rval == MB_SUCCESS);
339 if (MB_SUCCESS != rval) return rval;
340
341 //CartVect P0(&coords[0]);
342 //CartVect P3(&coords[3]);
343
344 //double normalVec[6];
345 CartVect N[2];
346 //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
347 rval = _mb->tag_get_data(_gradientTag, conn2, 2, (double*) &N[0]);
348 assert(rval == MB_SUCCESS);
349 if (MB_SUCCESS != rval) return rval;
350
351 CartVect T[2]; // T0, T3
352
353 rval = _mb->tag_get_data(_tangentsTag, &edge, 1, &T[0]);
354 assert(rval == MB_SUCCESS);
355 if (MB_SUCCESS != rval) return rval;
356
357 rval = init_edge_control_points(P[0], P[1], N[0], N[1], T[0], T[1], ctrl_pts);
358 assert(rval == MB_SUCCESS);
359 if (MB_SUCCESS != rval) return rval;
360
361 rval = _mb->tag_set_data(_edgeCtrlTag, &edge, 1, &ctrl_pts[0]);
362 assert(rval == MB_SUCCESS);
363 if (MB_SUCCESS != rval) return rval;
364
365 if (debug_surf_eval1)
366 {
367 std::cout << "edge: " << _mb-> id_from_handle(edge) << " tangents: "
368 << T[0] << T[1] << std::endl;
369 std::cout << " points: " << P[0] << " " << P[1] << std::endl;
370 std::cout << " normals: " << N[0] << " " << N[1] << std::endl;
371 std::cout << " Control points " << ctrl_pts[0] << " " << ctrl_pts[1]
372 << " " << ctrl_pts[2] << std::endl;
373 }
374
375 return MB_SUCCESS;
376 }
377
compute_tangents_for_each_edge()378 ErrorCode SmoothFace::compute_tangents_for_each_edge()
379 // they will be used for control points
380 {
381 double defTangents[6] = { 0., 0., 0., 0., 0., 0. };
382 ErrorCode rval = _mb->tag_get_handle("TANGENTS", 6, MB_TYPE_DOUBLE,
383 _tangentsTag, MB_TAG_DENSE|MB_TAG_CREAT, &defTangents);
384 if (MB_SUCCESS != rval)
385 return MB_FAILURE;
386
387 // now, compute Tangents for all edges that are not on boundary, so they are not marked
388 for (Range::iterator it = _edges.begin(); it != _edges.end(); ++it)
389 {
390 EntityHandle edg = *it;
391
392 int nnodes;
393 const EntityHandle * conn2;//
394 _mb->get_connectivity(edg, conn2, nnodes);
395 assert(nnodes == 2);
396 CartVect P[2]; // store the coordinates for the nodes
397 rval = _mb->get_coords(conn2, 2, (double *) &P[0]);
398 if (MB_SUCCESS != rval) return rval;
399 assert(rval==MB_SUCCESS);
400 CartVect T[2];
401 T[0] = P[1] - P[0];
402 T[0].normalize();
403 T[1] = T[0]; //
404 _mb->tag_set_data(_tangentsTag, &edg, 1, (double*) &T[0]);// set the tangents computed at every edge
405 }
406 return MB_SUCCESS;
407 }
408 // start copy
409 //===========================================================================
410 //Function Name: init_edge_control_points
411 //
412 //Member Type: PRIVATE
413 //Description: compute the control points for an edge
414 //===========================================================================
init_edge_control_points(CartVect & P0,CartVect & P3,CartVect & N0,CartVect & N3,CartVect & T0,CartVect & T3,CartVect * Pi)415 ErrorCode SmoothFace::init_edge_control_points(CartVect &P0, CartVect &P3,
416 CartVect &N0, CartVect &N3, CartVect &T0, CartVect &T3, CartVect * Pi)
417 {
418 CartVect Vi[4];
419 Vi[0] = P0;
420 Vi[3] = P3;
421 CartVect P03(P3 - P0);
422 double di = P03.length();
423 double ai = N0 % N3;// this is the dot operator, the same as in cgm for CubitVector
424 double ai0 = N0 % T0;
425 double ai3 = N3 % T3;
426 double denom = 4 - ai * ai;
427 if (fabs(denom) < 1e-20)
428 {
429 return MB_FAILURE; // CUBIT_FAILURE;
430 }
431 double row = 6.0e0 * (2.0e0 * ai0 + ai * ai3) / denom;
432 double omega = 6.0e0 * (2.0e0 * ai3 + ai * ai0) / denom;
433 Vi[1] = Vi[0] + (di * (((6.0e0 * T0) - ((2.0e0 * row) * N0) + (omega * N3))
434 / 18.0e0));
435 Vi[2] = Vi[3] - (di * (((6.0e0 * T3) + (row * N0) - ((2.0e0 * omega) * N3))
436 / 18.0e0));
437 //CartVect Wi[3];
438 //Wi[0] = Vi[1] - Vi[0];
439 //Wi[1] = Vi[2] - Vi[1];
440 //Wi[2] = Vi[3] - Vi[2];
441
442 Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
443 Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
444 Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
445
446 return MB_SUCCESS;
447 }
448
find_edges_orientations(EntityHandle edges[3],const EntityHandle * conn3,int orient[3])449 ErrorCode SmoothFace::find_edges_orientations(EntityHandle edges[3],
450 const EntityHandle * conn3, int orient[3])// maybe we will set it?
451 {
452 // find the edge that is adjacent to 2 vertices at a time
453 for (int i = 0; i < 3; i++)
454 {
455 // edge 0 is 1-2, 1 is 3-1, 2 is 0-1
456 EntityHandle v[2];
457 v[0] = conn3[(i + 1) % 3];
458 v[1] = conn3[(i + 2) % 3];
459 std::vector<EntityHandle> adjacencies;
460 // generate all edges for these two hexes
461 ErrorCode rval = _mb->get_adjacencies(v, 2, 1, false, adjacencies,
462 Interface::INTERSECT);
463 if (MB_SUCCESS != rval) return rval;
464
465 // find the edge connected to both vertices, and then see its orientation
466 assert(adjacencies.size() == 1);
467 const EntityHandle * conn2 = NULL;
468 int nnodes = 0;
469 rval = _mb->get_connectivity(adjacencies[0], conn2, nnodes);
470 assert(rval == MB_SUCCESS);
471 assert(2 == nnodes);
472 edges[i] = adjacencies[0];
473 // what is the story morning glory?
474 if (conn2[0] == v[0] && conn2[1] == v[1])
475 orient[i] = 1;
476 else if (conn2[0] == v[1] && conn2[1] == v[0])
477 orient[i] = -1;
478 else
479 return MB_FAILURE;
480 }
481 return MB_SUCCESS;
482 }
compute_internal_control_points_on_facets(double,Tag facetCtrlTag,Tag facetEdgeCtrlTag)483 ErrorCode SmoothFace::compute_internal_control_points_on_facets(double ,
484 Tag facetCtrlTag, Tag facetEdgeCtrlTag)
485 {
486 // collect from each triangle the control points in order
487 //
488
489 _facetCtrlTag = facetCtrlTag;
490 _facetEdgeCtrlTag = facetEdgeCtrlTag;
491
492 for (Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it)
493 {
494 EntityHandle tri = *it;
495 // first get connectivity, and the edges
496 // we need a fast method to retrieve the adjacent edges to each triangle
497 const EntityHandle * conn3;
498 int nnodes;
499 ErrorCode rval = _mb->get_connectivity(tri, conn3, nnodes);
500 assert(MB_SUCCESS == rval);
501 if (MB_SUCCESS != rval) return rval;
502 assert(3 == nnodes);
503
504 // would it be easier to do
505 CartVect vNode[3];// position at nodes
506 rval = _mb->get_coords(conn3, 3, (double*) &vNode[0]);
507 assert(MB_SUCCESS == rval);
508 if (MB_SUCCESS != rval) return rval;
509
510 // get gradients (normal) at each node of triangle
511 CartVect NN[3];
512 rval = _mb->tag_get_data(_gradientTag, conn3, 3, &NN[0]);
513 assert(MB_SUCCESS == rval);
514 if (MB_SUCCESS != rval) return rval;
515
516 EntityHandle edges[3];
517 int orient[3]; // + 1 or -1, if the edge is positive or negative within the face
518 rval = find_edges_orientations(edges, conn3, orient);// maybe we will set it?
519 assert(MB_SUCCESS == rval);
520 if (MB_SUCCESS != rval) return rval;
521 // maybe we will store some tags with edges and their orientation with respect to
522 // a triangle;
523 CartVect P[3][5];
524 CartVect N[6], G[6];
525 // create the linear array for control points on edges, for storage (expensive !!!)
526 CartVect CP[9];
527 int index = 0;
528 // maybe store a tag / entity handle for edges?
529 for (int i = 0; i < 3; i++)
530 {
531 // populate P and N with the right vectors
532 int i1 = (i + 1) % 3; // the first node of the edge
533 int i2 = (i + 2) % 3; // the second node of the edge
534 N[2 * i] = NN[i1];
535 N[2 * i + 1] = NN[i2];
536 P[i][0] = vNode[i1];
537 rval = _mb->tag_get_data(_edgeCtrlTag, &edges[i], 1, &(P[i][1]));
538 // if sense is -1, swap 1 and 3 control points
539 if (orient[i] == -1)
540 {
541 CartVect tmp;
542 tmp = P[i][1];
543 P[i][1] = P[i][3];
544 P[i][3] = tmp;
545 }
546 P[i][4] = vNode[i2];
547 for (int j = 1; j < 4; j++)
548 CP[index++] = P[i][j];
549
550 // the first edge control points
551 }
552
553 // stat = facet->get_edge_control_points( P );
554 init_facet_control_points(N, P, G);
555 // what do we need to store in the tag control points?
556 rval = _mb->tag_set_data(_facetCtrlTag, &tri, 1, &G[0]);
557 assert(MB_SUCCESS == rval);
558 if (MB_SUCCESS != rval) return rval;
559
560 // store here again the 9 control points on the edges
561 rval = _mb->tag_set_data(_facetEdgeCtrlTag, &tri, 1, &CP[0]);
562 assert(MB_SUCCESS == rval);
563 if (MB_SUCCESS != rval) return rval;
564 // look at what we retrieve later
565
566 // adjust the bounding box
567 int j = 0;
568 for (j = 0; j < 3; j++)
569 adjust_bounding_box(vNode[j]);
570 // edge control points
571 for (j = 0; j < 9; j++)
572 adjust_bounding_box(CP[j]);
573 // internal facet control points
574 for (j = 0; j < 6; j++)
575 adjust_bounding_box(G[j]);
576
577 }
578 return MB_SUCCESS;
579 }
adjust_bounding_box(CartVect & vect)580 void SmoothFace::adjust_bounding_box(CartVect & vect)
581 {
582 // _minim, _maxim
583 for (int j = 0; j < 3; j++)
584 {
585 if (_minim[j] > vect[j])
586 _minim[j] = vect[j];
587 if (_maxim[j] < vect[j])
588 _maxim[j] = vect[j];
589 }
590 }
591 //===============================================================
592 ////Function Name: init_facet_control_points
593 ////
594 ////Member Type: PRIVATE
595 ////Description: compute the control points for a facet
596 ////===============================================================
init_facet_control_points(CartVect N[6],CartVect P[3][5],CartVect G[6])597 ErrorCode SmoothFace::init_facet_control_points(CartVect N[6], // vertex normals (per edge)
598 CartVect P[3][5], // edge control points
599 CartVect G[6]) // return internal control points
600 {
601 CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
602 double denom;
603 double lambda[2], mu[2];
604
605 ErrorCode rval = MB_SUCCESS;
606
607 for (int i = 0; i < 3; i++)
608 {
609 N0 = N[i * 2];
610 N3 = N[i * 2 + 1];
611 Vi[0] = P[i][0];
612 Vi[1] = (P[i][1] - 0.25 * P[i][0]) / 0.75;
613 Vi[2] = (P[i][3] - 0.25 * P[i][4]) / 0.75;
614 Vi[3] = P[i][4];
615 Wi[0] = Vi[1] - Vi[0];
616 Wi[1] = Vi[2] - Vi[1];
617 Wi[2] = Vi[3] - Vi[2];
618 Di[0] = P[(i + 2) % 3][3] - 0.5 * (P[i][1] + P[i][0]);
619 Di[3] = P[(i + 1) % 3][1] - 0.5 * (P[i][4] + P[i][3]);
620 Ai[0] = (N0 * Wi[0]) / Wi[0].length();
621 Ai[2] = (N3 * Wi[2]) / Wi[2].length();
622 Ai[1] = Ai[0] + Ai[2];
623 denom = Ai[1].length();
624 Ai[1] /= denom;
625 lambda[0] = (Di[0] % Wi[0]) / (Wi[0] % Wi[0]);
626 lambda[1] = (Di[3] % Wi[2]) / (Wi[2] % Wi[2]);
627 mu[0] = (Di[0] % Ai[0]);
628 mu[1] = (Di[3] % Ai[2]);
629 G[i * 2] = 0.5 * (P[i][1] + P[i][2]) + 0.66666666666666 * lambda[0] * Wi[1]
630 + 0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0]
631 * Ai[1] + 0.33333333333333 * mu[1] * Ai[0];
632 G[i * 2 + 1] = 0.5 * (P[i][2] + P[i][3]) + 0.33333333333333 * lambda[0]
633 * Wi[2] + 0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333
634 * mu[0] * Ai[2] + 0.66666666666666 * mu[1] * Ai[1];
635 }
636 return rval;
637 }
638
DumpModelControlPoints()639 void SmoothFace::DumpModelControlPoints()
640 {
641 // here, we will dump all control points from edges and facets (6 control points for each facet)
642 // we may also create some edges; maybe later...
643 // create a point3D file
644 // output a Point3D file (special visit format)
645 unsigned long setId = _mb->id_from_handle(_set);
646 char name[50] = { 0 };
647 sprintf(name, "%lucontrol.Point3D", setId);// name should be something 2control.Point3D
648 std::ofstream point3DFile;
649 point3DFile.open(name);//("control.Point3D");
650 point3DFile << "# x y z \n";
651 std::ofstream point3DEdgeFile;
652 sprintf(name, "%lucontrolEdge.Point3D", setId);//
653 point3DEdgeFile.open(name);//("controlEdge.Point3D");
654 point3DEdgeFile << "# x y z \n";
655 std::ofstream smoothPoints;
656 sprintf(name, "%lusmooth.Point3D", setId);//
657 smoothPoints.open(name);//("smooth.Point3D");
658 smoothPoints << "# x y z \n";
659 CartVect controlPoints[3]; // edge control points
660 for (Range::iterator it = _edges.begin(); it != _edges.end(); ++it)
661 {
662 EntityHandle edge = *it;
663
664 _mb->tag_get_data(_edgeCtrlTag, &edge, 1, (double*) &controlPoints[0]);
665 for (int i = 0; i < 3; i++)
666 {
667 CartVect & c = controlPoints[i];
668 point3DEdgeFile << std::setprecision(11) << c[0] << " " << c[1] << " "
669 << c[2] << " \n";
670 }
671 }
672 CartVect controlTriPoints[6]; // triangle control points
673 CartVect P_facet[3];// result in 3 "mid" control points
674 for (Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); ++it2)
675 {
676 EntityHandle tri = *it2;
677
678 _mb->tag_get_data(_facetCtrlTag, &tri, 1, (double*) &controlTriPoints[0]);
679
680 // draw a line of points between pairs of control points
681 int numPoints = 7;
682 for (int n = 0; n < numPoints; n++)
683 {
684 double a = 1. * n / (numPoints - 1);
685 double b = 1.0 - a;
686
687 P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
688 //1,2,1
689 P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
690 //1,1,2
691 P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
692 for (int i = 0; i < 3; i++)
693 {
694 CartVect & c = P_facet[i];
695 point3DFile << std::setprecision(11) << c[0] << " " << c[1] << " "
696 << c[2] << " \n";
697 }
698 }
699
700 // evaluate for each triangle a lattice of points
701 int N = 40;
702 for (int k = 0; k <= N; k++)
703 {
704 for (int m = 0; m <= N - k; m++)
705 {
706 int n = N - m - k;
707 CartVect areacoord(1. * k / N, 1. * m / N, 1. * n / N);
708 CartVect pt;
709 eval_bezier_patch(tri, areacoord, pt);
710 smoothPoints << std::setprecision(11) << pt[0] << " " << pt[1] << " "
711 << pt[2] << " \n";
712
713 }
714 }
715 }
716 point3DFile.close();
717 smoothPoints.close();
718 point3DEdgeFile.close();
719 return;
720 }
721 //===========================================================================
722 //Function Name: evaluate_single
723 //
724 //Member Type: PUBLIC
725 //Description: evaluate edge not associated with a facet (this is used
726 // by camal edge mesher!!!)
727 //Note: t is a value from 0 to 1, for us
728 //===========================================================================
evaluate_smooth_edge(EntityHandle eh,double & tt,CartVect & outv)729 ErrorCode SmoothFace::evaluate_smooth_edge(EntityHandle eh, double &tt,
730 CartVect & outv)
731 {
732 CartVect P[2]; // P0 and P1
733 CartVect controlPoints[3]; // edge control points
734 double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
735
736 // project the position to the linear edge
737 // t is from 0 to 1 only!!
738 //double tt = (t + 1) * 0.5;
739 if (tt <= 0.0)
740 tt = 0.0;
741 if (tt >= 1.0)
742 tt = 1.0;
743
744 int nnodes = 0;
745 const EntityHandle * conn2 = NULL;
746 ErrorCode rval = _mb->get_connectivity(eh, conn2, nnodes);
747 assert(rval == MB_SUCCESS);
748 if (MB_SUCCESS != rval) return rval;
749
750 rval = _mb->get_coords(conn2, 2, (double*) &P[0]);
751 assert(rval == MB_SUCCESS);
752 if (MB_SUCCESS != rval) return rval;
753
754 rval = _mb->tag_get_data(_edgeCtrlTag, &eh, 1, (double*) &controlPoints[0]);
755 assert(rval == MB_SUCCESS);
756 if (MB_SUCCESS != rval) return rval;
757
758 t2 = tt * tt;
759 t3 = t2 * tt;
760 t4 = t3 * tt;
761 one_minus_t = 1. - tt;
762 one_minus_t2 = one_minus_t * one_minus_t;
763 one_minus_t3 = one_minus_t2 * one_minus_t;
764 one_minus_t4 = one_minus_t3 * one_minus_t;
765
766 outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6.
767 * one_minus_t2 * t2 * controlPoints[1] + 4. * one_minus_t * t3
768 * controlPoints[2] + t4 * P[1];
769
770 return MB_SUCCESS;
771 }
772
eval_bezier_patch(EntityHandle tri,CartVect & areacoord,CartVect & pt)773 ErrorCode SmoothFace::eval_bezier_patch(EntityHandle tri, CartVect &areacoord,
774 CartVect &pt)
775 {
776 //
777 // interpolate internal control points
778
779 CartVect gctrl_pts[6];
780 // get the control points facet->get_control_points( gctrl_pts );
781 //init_facet_control_points( N, P, G) ;
782 // what do we need to store in the tag control points?
783 ErrorCode rval = _mb->tag_get_data(_facetCtrlTag, &tri, 1, &gctrl_pts[0]);// get all 6 control points
784 assert(MB_SUCCESS == rval);
785 if (MB_SUCCESS != rval) return rval;
786 const EntityHandle * conn3 = NULL;
787 int nnodes = 0;
788 rval = _mb->get_connectivity(tri, conn3, nnodes);
789 assert(MB_SUCCESS == rval);
790
791 CartVect vN[3];
792 _mb->get_coords(conn3, 3, (double*) &vN[0]); // fill the coordinates of the vertices
793
794 if (fabs(areacoord[1] + areacoord[2]) < 1.0e-6)
795 {
796 pt = vN[0];
797 return MB_SUCCESS;
798 }
799 if (fabs(areacoord[0] + areacoord[2]) < 1.0e-6)
800 {
801 pt = vN[0];
802 return MB_SUCCESS;
803 }
804 if (fabs(areacoord[0] + areacoord[1]) < 1.0e-6)
805 {
806 pt = vN[0];
807 return MB_SUCCESS;
808 }
809
810 CartVect P_facet[3];
811 //2,1,1
812 P_facet[0] = (1.0e0 / (areacoord[1] + areacoord[2])) * (areacoord[1]
813 * gctrl_pts[3] + areacoord[2] * gctrl_pts[4]);
814 //1,2,1
815 P_facet[1] = (1.0e0 / (areacoord[0] + areacoord[2])) * (areacoord[0]
816 * gctrl_pts[0] + areacoord[2] * gctrl_pts[5]);
817 //1,1,2
818 P_facet[2] = (1.0e0 / (areacoord[0] + areacoord[1])) * (areacoord[0]
819 * gctrl_pts[1] + areacoord[1] * gctrl_pts[2]);
820
821 // sum the contribution from each of the control points
822
823 pt = CartVect(0.);// set all to 0, we start adding / accumulating different parts
824 // first edge is from node 0 to 1, index 2 in
825
826 // retrieve the points, in order, and the control points on edges
827
828 // store here again the 9 control points on the edges
829 CartVect CP[9];
830 rval = _mb->tag_get_data(_facetEdgeCtrlTag, &tri, 1, &CP[0]);
831 assert(MB_SUCCESS == rval);
832
833 //CubitFacetEdge *edge;
834 //edge = facet->edge(2);! start with edge 2, from 0-1
835 int k = 0;
836 CartVect ctrl_pts[5];
837 //edge->control_points(facet, ctrl_pts);
838 ctrl_pts[0] = vN[0]; //
839 for (k = 1; k < 4; k++)
840 ctrl_pts[k] = CP[k + 5]; // for edge index 2
841 ctrl_pts[4] = vN[1]; //
842
843 //i=4; j=0; k=0;
844 double B = mbquart(areacoord[0]);
845 pt += B * ctrl_pts[0];
846
847 //i=3; j=1; k=0;
848 B = 4.0 * mbcube(areacoord[0]) * areacoord[1];
849 pt += B * ctrl_pts[1];
850
851 //i=2; j=2; k=0;
852 B = 6.0 * mbsqr(areacoord[0]) * mbsqr(areacoord[1]);
853 pt += B * ctrl_pts[2];
854
855 //i=1; j=3; k=0;
856 B = 4.0 * areacoord[0] * mbcube(areacoord[1]);
857 pt += B * ctrl_pts[3];
858
859 //edge = facet->edge(0);
860 //edge->control_points(facet, ctrl_pts);
861 // edge index 0, from 1 to 2
862 ctrl_pts[0] = vN[1]; //
863 for (k = 1; k < 4; k++)
864 ctrl_pts[k] = CP[k - 1]; // for edge index 0
865 ctrl_pts[4] = vN[2]; //
866
867 //i=0; j=4; k=0;
868 B = mbquart(areacoord[1]);
869 pt += B * ctrl_pts[0];
870
871 //i=0; j=3; k=1;
872 B = 4.0 * mbcube(areacoord[1]) * areacoord[2];
873 pt += B * ctrl_pts[1];
874
875 //i=0; j=2; k=2;
876 B = 6.0 * mbsqr(areacoord[1]) * mbsqr(areacoord[2]);
877 pt += B * ctrl_pts[2];
878
879 //i=0; j=1; k=3;
880 B = 4.0 * areacoord[1] * mbcube(areacoord[2]);
881 pt += B * ctrl_pts[3];
882
883 //edge = facet->edge(1);
884 //edge->control_points(facet, ctrl_pts);
885 // edge index 1, from 2 to 0
886 ctrl_pts[0] = vN[2]; //
887 for (k = 1; k < 4; k++)
888 ctrl_pts[k] = CP[k + 2]; // for edge index 0
889 ctrl_pts[4] = vN[0]; //
890
891 //i=0; j=0; k=4;
892 B = mbquart(areacoord[2]);
893 pt += B * ctrl_pts[0];
894
895 //i=1; j=0; k=3;
896 B = 4.0 * areacoord[0] * mbcube(areacoord[2]);
897 pt += B * ctrl_pts[1];
898
899 //i=2; j=0; k=2;
900 B = 6.0 * mbsqr(areacoord[0]) * mbsqr(areacoord[2]);
901 pt += B * ctrl_pts[2];
902
903 //i=3; j=0; k=1;
904 B = 4.0 * mbcube(areacoord[0]) * areacoord[2];
905 pt += B * ctrl_pts[3];
906
907 //i=2; j=1; k=1;
908 B = 12.0 * mbsqr(areacoord[0]) * areacoord[1] * areacoord[2];
909 pt += B * P_facet[0];
910
911 //i=1; j=2; k=1;
912 B = 12.0 * areacoord[0] * mbsqr(areacoord[1]) * areacoord[2];
913 pt += B * P_facet[1];
914
915 //i=1; j=1; k=2;
916 B = 12.0 * areacoord[0] * areacoord[1] * mbsqr(areacoord[2]);
917 pt += B * P_facet[2];
918
919 return MB_SUCCESS;
920
921 }
922
923 //===========================================================================
924 //Function Name: project_to_facet_plane
925 //
926 //Member Type: PUBLIC
927 //Descriptoin: Project a point to the plane of a facet
928 //===========================================================================
project_to_facet_plane(EntityHandle tri,CartVect & pt,CartVect & point_on_plane,double & dist_to_plane)929 void SmoothFace::project_to_facet_plane(EntityHandle tri, CartVect &pt,
930 CartVect &point_on_plane, double &dist_to_plane)
931 {
932 double plane[4];
933
934 ErrorCode rval = _mb->tag_get_data(_planeTag, &tri, 1, plane);
935 if (MB_SUCCESS != rval) return;
936 assert(rval == MB_SUCCESS);
937 // _planeTag
938 CartVect normal(&plane[0]);// just first 3 components are used
939
940 double dist = normal % pt + plane[3]; // coeff d is saved!!!
941 dist_to_plane = fabs(dist);
942 point_on_plane = pt - dist * normal;
943
944 return;
945 }
946
947 //===========================================================================
948 //Function Name: facet_area_coordinate
949 //
950 //Member Type: PUBLIC
951 //Descriptoin: Determine the area coordinates of a point on the plane
952 // of a facet
953 //===========================================================================
facet_area_coordinate(EntityHandle facet,CartVect & pt_on_plane,CartVect & areacoord)954 void SmoothFace::facet_area_coordinate(EntityHandle facet,
955 CartVect & pt_on_plane, CartVect & areacoord)
956 {
957
958 const EntityHandle * conn3 = NULL;
959 int nnodes = 0;
960 ErrorCode rval = _mb->get_connectivity(facet, conn3, nnodes);
961 assert(MB_SUCCESS == rval);
962 if (rval) {} // empty statement to prevent compiler warning
963
964 //double coords[9]; // store the coordinates for the nodes
965 //_mb->get_coords(conn3, 3, coords);
966 CartVect p[3];
967 rval = _mb->get_coords(conn3, 3, (double*) &p[0]);
968 assert(MB_SUCCESS == rval);
969 if (rval) {} // empty statement to prevent compiler warning
970 double plane[4];
971
972 rval = _mb->tag_get_data(_planeTag, &facet, 1, plane);
973 assert(rval == MB_SUCCESS);
974 if (rval) {} // empty statement to prevent compiler warning
975 CartVect normal(&plane[0]);// just first 3 components are used
976
977 double area2;
978
979 double tol = GEOMETRY_RESABS * 1.e-5;// 1.e-11;
980
981 CartVect v1(p[1] - p[0]);
982 CartVect v2(p[2] - p[0]);
983
984 area2 = (v1 * v2).length_squared();// the same for CartVect
985 if (area2 < 100 * tol)
986 {
987 tol = .01 * area2;
988 }
989 CartVect absnorm(fabs(normal[0]), fabs(normal[1]), fabs(normal[2]));
990
991 // project to the closest coordinate plane so we only have to do this in 2D
992
993 if (absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2])
994 {
995 area2 = determ3( p[0][1], p[0][2],
996 p[1][1], p[1][2],
997 p[2][1], p[2][2]);
998 if (fabs(area2) < tol)
999 {
1000 areacoord = CartVect(-std::numeric_limits<double>::min());// .set( -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min() );
1001 }
1002 else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
1003 {
1004 areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
1005 }
1006 else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
1007 {
1008 areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
1009 }
1010 else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
1011 {
1012 areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
1013 }
1014 else
1015 {
1016
1017 areacoord[0] = determ3(pt_on_plane[1], pt_on_plane[2],
1018 p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
1019
1020 areacoord[1] = determ3( p[0][1], p[0][2],
1021 pt_on_plane[1], pt_on_plane[2],
1022 p[2][1], p[2][2]) / area2;
1023
1024 areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2],
1025 pt_on_plane[1], pt_on_plane[2]) / area2;
1026 }
1027 }
1028 else if (absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2])
1029 {
1030
1031 area2 = determ3(p[0][0], p[0][2],
1032 p[1][0], p[1][2],
1033 p[2][0], p[2][2]);
1034 if (fabs(area2) < tol)
1035 {
1036 areacoord = CartVect(-std::numeric_limits<double>::min());//.set( -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min() );
1037 }
1038 else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
1039 {
1040 areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
1041 }
1042 else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
1043 {
1044 areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
1045 }
1046 else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
1047 {
1048 areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
1049 }
1050 else
1051 {
1052
1053 areacoord[0] = determ3(pt_on_plane[0], pt_on_plane[2],
1054 p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
1055
1056 areacoord[1] = determ3( p[0][0], p[0][2],
1057 pt_on_plane[0], pt_on_plane[2],
1058 p[2][0], p[2][2]) / area2;
1059
1060 areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2],
1061 pt_on_plane[0], pt_on_plane[2]) / area2;
1062
1063 }
1064 }
1065 else
1066 {
1067 /*area2 = determ3(pt0->x(), pt0->y(),
1068 pt1->x(), pt1->y(),
1069 pt2->x(), pt2->y());*/
1070 area2 = determ3( p[0][0], p[0][1],
1071 p[1][0], p[1][1],
1072 p[2][0], p[2][1]);
1073 if (fabs(area2) < tol)
1074 {
1075 areacoord = CartVect(-std::numeric_limits<double>::min());//.set( -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min(), -std::numeric_limits<double>::min() );
1076 }
1077 else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
1078 {
1079 areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
1080 }
1081 else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
1082 {
1083 areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
1084 }
1085 else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
1086 {
1087 areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
1088 }
1089 else
1090 {
1091
1092 areacoord[0] = determ3(pt_on_plane[0], pt_on_plane[1],
1093 p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
1094
1095 areacoord[1] = determ3( p[0][0], p[0][1],
1096 pt_on_plane[0], pt_on_plane[1],
1097 p[2][0], p[2][1]) / area2;
1098
1099 areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1],
1100 pt_on_plane[0], pt_on_plane[1]) / area2;
1101 }
1102 }
1103 }
1104
project_to_facets_main(CartVect & this_point,bool trim,bool & outside,CartVect * closest_point_ptr,CartVect * normal_ptr)1105 ErrorCode SmoothFace::project_to_facets_main(CartVect &this_point, bool trim,
1106 bool & outside, CartVect * closest_point_ptr, CartVect * normal_ptr)
1107 {
1108
1109 // if there are a lot of facets on this surface - use the OBB search first
1110 // to narrow the selection
1111
1112 _evaluationsCounter++;
1113 double tolerance = 1.e-5;
1114 std::vector<EntityHandle> facets_out;
1115 // we will start with a list of facets anyway, the best among them wins
1116 ErrorCode rval = _my_geomTopoTool->obb_tree()->closest_to_location(
1117 (double*) &this_point, _obb_root, tolerance, facets_out);
1118 if (MB_SUCCESS != rval) return rval;
1119
1120 int interpOrder = 4;
1121 double compareTol = 1.e-5;
1122 EntityHandle lastFacet = facets_out.front();
1123 rval = project_to_facets(facets_out, lastFacet, interpOrder, compareTol,
1124 this_point, trim, outside, closest_point_ptr, normal_ptr);
1125
1126 return rval;
1127 }
project_to_facets(std::vector<EntityHandle> & facet_list,EntityHandle & lastFacet,int interpOrder,double compareTol,CartVect & this_point,bool,bool & outside,CartVect * closest_point_ptr,CartVect * normal_ptr)1128 ErrorCode SmoothFace::project_to_facets(std::vector<EntityHandle> & facet_list,
1129 EntityHandle & lastFacet, int interpOrder, double compareTol,
1130 CartVect &this_point, bool , bool & outside,
1131 CartVect *closest_point_ptr, CartVect * normal_ptr)
1132 {
1133
1134 bool outside_facet = false;
1135 bool best_outside_facet = true;
1136 double mindist = 1.e20;
1137 CartVect close_point, best_point(mindist, mindist, mindist), best_areacoord;
1138 EntityHandle best_facet = 0L;// no best facet found yet
1139 EntityHandle facet;
1140 assert(facet_list.size() > 0);
1141
1142 double big_dist = compareTol * 1.0e3;
1143
1144 // from the list of close facets, determine the closest point
1145 for (size_t i = 0; i < facet_list.size(); i++)
1146 {
1147 facet = facet_list[i];
1148 CartVect pt_on_plane;
1149 double dist_to_plane;
1150 project_to_facet_plane(facet, this_point, pt_on_plane, dist_to_plane);
1151
1152 CartVect areacoord;
1153 //CartVect close_point;
1154 facet_area_coordinate(facet, pt_on_plane, areacoord);
1155 if (interpOrder != 0)
1156 {
1157
1158 // modify the areacoord - project to the bezier patch- snaps to the
1159 // edge of the patch if necessary
1160
1161
1162 if (project_to_facet(facet, this_point, areacoord, close_point,
1163 outside_facet, compareTol) != MB_SUCCESS)
1164 {
1165 return MB_FAILURE;
1166 }
1167 //if (closest_point_ptr)
1168 //*closest_point_ptr = close_point;
1169 }
1170 // keep track of the minimum distance
1171
1172 double dist = (close_point - this_point).length();//close_point.distance_between(this_point);
1173 if ((best_outside_facet == outside_facet && dist < mindist)
1174 || (best_outside_facet && !outside_facet && (dist < big_dist
1175 || best_facet == 0L/*!best_facet*/)))
1176 {
1177 mindist = dist;
1178 best_point = close_point;
1179 best_facet = facet;
1180 best_areacoord = areacoord;
1181 best_outside_facet = outside_facet;
1182
1183 if (dist < compareTol)
1184 {
1185 break;
1186 }
1187 big_dist = 10.0 * mindist;
1188 }
1189 //facet->marked(1);
1190 //used_facet_list.append(facet);
1191
1192 }
1193
1194 if (normal_ptr)
1195 {
1196 CartVect normal;
1197 if (eval_bezier_patch_normal(best_facet, best_areacoord, normal)
1198 != MB_SUCCESS)
1199 {
1200 return MB_FAILURE;
1201 }
1202 *normal_ptr = normal;
1203 }
1204
1205 if (closest_point_ptr)
1206 {
1207 *closest_point_ptr = best_point;
1208 }
1209
1210 outside = best_outside_facet;
1211 lastFacet = best_facet;
1212
1213 return MB_SUCCESS;
1214 //end copy
1215 }
1216
1217 //===========================================================================
1218 //Function Name: project_to_patch
1219 //
1220 //Member Type: PUBLIC
1221 //Description: Project a point to a bezier patch. Pass in the areacoord
1222 // of the point projected to the linear facet. Function
1223 // assumes that the point is contained within the patch -
1224 // if not, it will project to one of its edges.
1225 //===========================================================================
project_to_patch(EntityHandle facet,CartVect & ac,CartVect & pt,CartVect & eval_pt,CartVect * eval_norm,bool & outside,double compare_tol,int edge_id)1226 ErrorCode SmoothFace::project_to_patch(EntityHandle facet, // (IN) the facet where the patch is defined
1227 CartVect &ac, // (IN) area coordinate initial guess (from linear facet)
1228 CartVect &pt, // (IN) point we are projecting to patch
1229 CartVect &eval_pt, // (OUT) The projected point
1230 CartVect *eval_norm, // (OUT) normal at evaluated point
1231 bool &outside, // (OUT) the closest point on patch to pt is on an edge
1232 double compare_tol, // (IN) comparison tolerance
1233 int edge_id) // (IN) only used if this is to be projected to one
1234 // of the edges. Otherwise, should be -1
1235 {
1236 ErrorCode status = MB_SUCCESS;
1237
1238 // see if we are at a vertex
1239
1240 #define INCR 0.01
1241 const double tol = compare_tol;
1242
1243 if (is_at_vertex(facet, pt, ac, compare_tol, eval_pt, eval_norm))
1244 {
1245 outside = false;
1246 return MB_SUCCESS;
1247 }
1248
1249 // check if the start ac is inside the patch -if not, then move it there
1250
1251 int nout = 0;
1252 const double atol = 0.001;
1253 if (move_ac_inside(ac, atol))
1254 nout++;
1255
1256 int diverge = 0;
1257 int iter = 0;
1258 CartVect newpt;
1259 eval_bezier_patch(facet, ac, newpt);
1260 CartVect move = pt - newpt;
1261 double lastdist = move.length();
1262 double bestdist = lastdist;
1263 CartVect bestac = ac;
1264 CartVect bestpt = newpt;
1265 CartVect bestnorm(0, 0, 0);
1266
1267 // If we are already close enough, then return now
1268
1269 if (lastdist <= tol && !eval_norm && nout == 0)
1270 {
1271 eval_pt = pt;
1272 outside = false;
1273 return status;
1274 }
1275
1276 double ratio, mag, umove, vmove, det, distnew, movedist;
1277 CartVect lastpt = newpt;
1278 CartVect lastac = ac;
1279 CartVect norm;
1280 CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
1281 CartVect du, dv, newac;
1282 bool done = false;
1283 while (!done)
1284 {
1285
1286 // We will be locating the projected point within the u,v,w coordinate
1287 // system of the triangular bezier patch. Since u+v+w=1, the components
1288 // are not linearly independent. We will choose only two of the
1289 // coordinates to use and compute the third.
1290
1291 int system;
1292 if (lastac[0] >= lastac[1] && lastac[0] >= lastac[2])
1293 {
1294 system = 0;
1295 }
1296 else if (lastac[1] >= lastac[2])
1297 {
1298 system = 1;
1299 }
1300 else
1301 {
1302 system = 2;
1303 }
1304
1305 // compute the surface derivatives with respect to each
1306 // of the barycentric coordinates
1307
1308
1309 if (system == 1 || system == 2)
1310 {
1311 xac[0] = lastac[0] + INCR; // xac.x( lastac.x() + INCR );
1312 if (lastac[1] + lastac[2] == 0.0)
1313 return MB_FAILURE;
1314 ratio = lastac[2] / (lastac[1] + lastac[2]);//ratio = lastac.z() / (lastac.y() + lastac.z());
1315 xac[1] = (1.0 - xac[0]) * (1.0 - ratio);//xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
1316 xac[2] = 1.0 - xac[0] - xac[1]; // xac.z( 1.0 - xac.x() - xac.y() );
1317 eval_bezier_patch(facet, xac, xpt);
1318 xvec = xpt - lastpt;
1319 xvec /= INCR;
1320 }
1321 if (system == 0 || system == 2)
1322 {
1323 yac[1] = (lastac[1] + INCR);//yac.y( lastac.y() + INCR );
1324 if (lastac[0] + lastac[2] == 0.0)//if (lastac.x() + lastac.z() == 0.0)
1325 return MB_FAILURE;
1326 ratio = lastac[2] / (lastac[0] + lastac[2]);//ratio = lastac.z() / (lastac.x() + lastac.z());
1327 yac[0] = ((1.0 - yac[1]) * (1.0 - ratio));//yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
1328 yac[2] = (1.0 - yac[0] - yac[1]);//yac.z( 1.0 - yac.x() - yac.y() );
1329 eval_bezier_patch(facet, yac, ypt);
1330 yvec = ypt - lastpt;
1331 yvec /= INCR;
1332 }
1333 if (system == 0 || system == 1)
1334 {
1335 zac[2] = (lastac[2] + INCR);//zac.z( lastac.z() + INCR );
1336 if (lastac[0] + lastac[1] == 0.0)//if (lastac.x() + lastac.y() == 0.0)
1337 return MB_FAILURE;
1338 ratio = lastac[1] / (lastac[0] + lastac[1]);//ratio = lastac.y() / (lastac.x() + lastac.y());
1339 zac[0] = ((1.0 - zac[2]) * (1.0 - ratio));//zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
1340 zac[1] = (1.0 - zac[0] - zac[2]);//zac.y( 1.0 - zac.x() - zac.z() );
1341 eval_bezier_patch(facet, zac, zpt);
1342 zvec = zpt - lastpt;
1343 zvec /= INCR;
1344 }
1345
1346 // compute the surface normal
1347
1348 switch (system) {
1349 case 0:
1350 du = yvec;
1351 dv = zvec;
1352 break;
1353 case 1:
1354 du = zvec;
1355 dv = xvec;
1356 break;
1357 case 2:
1358 du = xvec;
1359 dv = yvec;
1360 break;
1361 }
1362 norm = du * dv;
1363 mag = norm.length();
1364 if (mag < DBL_EPSILON)
1365 {
1366 return MB_FAILURE;
1367 // do something else here (it is likely a flat triangle -
1368 // so try evaluating just an edge of the bezier patch)
1369 }
1370 norm /= mag;
1371 if (iter == 0)
1372 bestnorm = norm;
1373
1374 // project the move vector to the tangent plane
1375
1376 move = (norm * move) * norm;
1377
1378 // compute an equivalent u-v-w vector
1379
1380 CartVect absnorm(fabs(norm[0]), fabs(norm[1]), fabs(norm[2]));
1381 if (absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0])
1382 {
1383 det = du[0] * dv[1] - dv[0] * du[1];
1384 if (fabs(det) <= DBL_EPSILON)
1385 {
1386 return MB_FAILURE; // do something else here
1387 }
1388 umove = (move[0] * dv[1] - dv[0] * move[1]) / det;
1389 vmove = (du[0] * move[1] - move[0] * du[1]) / det;
1390 }
1391 else if (absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0])
1392 {
1393 det = du[0] * dv[2] - dv[0] * du[2];
1394 if (fabs(det) <= DBL_EPSILON)
1395 {
1396 return MB_FAILURE;
1397 }
1398 umove = (move[0] * dv[2] - dv[0] * move[2]) / det;
1399 vmove = (du[0] * move[2] - move[0] * du[2]) / det;
1400 }
1401 else
1402 {
1403 det = du[1] * dv[2] - dv[1] * du[2];
1404 if (fabs(det) <= DBL_EPSILON)
1405 {
1406 return MB_FAILURE;
1407 }
1408 umove = (move[1] * dv[2] - dv[1] * move[2]) / det;
1409 vmove = (du[1] * move[2] - move[1] * du[2]) / det;
1410 }
1411
1412 /* === compute the new u-v coords and evaluate surface at new location */
1413
1414 switch (system) {
1415 case 0:
1416 newac[1] = (lastac[1] + umove);//newac.y( lastac.y() + umove );
1417 newac[2] = (lastac[2] + vmove);//newac.z( lastac.z() + vmove );
1418 newac[0] = (1.0 - newac[1] - newac[2]);//newac.x( 1.0 - newac.y() - newac.z() );
1419 break;
1420 case 1:
1421 newac[2] = (lastac[2] + umove);//newac.z( lastac.z() + umove );
1422 newac[0] = (lastac[0] + vmove);//newac.x( lastac.x() + vmove );
1423 newac[1] = (1.0 - newac[2] - newac[0]);//newac.y( 1.0 - newac.z() - newac.x() );
1424 break;
1425 case 2:
1426 newac[0] = (lastac[0] + umove);//newac.x( lastac.x() + umove );
1427 newac[1] = (lastac[1] + vmove);//newac.y( lastac.y() + vmove );
1428 newac[2] = (1.0 - newac[0] - newac[1]);//newac.z( 1.0 - newac.x() - newac.y() );
1429 break;
1430 }
1431
1432 // Keep it inside the patch
1433
1434 if (newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol)
1435 {
1436 nout = 0;
1437 }
1438 else
1439 {
1440 if (move_ac_inside(newac, atol))
1441 nout++;
1442 }
1443
1444 // Evaluate at the new location
1445
1446 if (edge_id != -1)
1447 ac_at_edge(newac, newac, edge_id); // move to edge first
1448 eval_bezier_patch(facet, newac, newpt);
1449
1450 // Check for convergence
1451
1452 distnew = (pt - newpt).length();//pt.distance_between(newpt);
1453 move = newpt - lastpt;
1454 movedist = move.length();
1455 if (movedist < tol || distnew < tol)
1456 {
1457 done = true;
1458 if (distnew < bestdist)
1459 {
1460 bestdist = distnew;
1461 bestac = newac;
1462 bestpt = newpt;
1463 bestnorm = norm;
1464 }
1465 }
1466 else
1467 {
1468
1469 // don't allow more than 30 iterations
1470
1471 iter++;
1472 if (iter > 30)
1473 {
1474 //if (movedist > tol * 100.0) nout=1;
1475 done = true;
1476 }
1477
1478 // Check for divergence - don't allow more than 5 divergent
1479 // iterations
1480
1481 if (distnew > lastdist)
1482 {
1483 diverge++;
1484 if (diverge > 10)
1485 {
1486 done = true;
1487 //if (movedist > tol * 100.0) nout=1;
1488 }
1489 }
1490
1491 // Check if we are continuing to project outside the facet.
1492 // If so, then stop now
1493
1494 if (nout > 3)
1495 {
1496 done = true;
1497 }
1498
1499 // set up for next iteration
1500
1501 if (!done)
1502 {
1503 if (distnew < bestdist)
1504 {
1505 bestdist = distnew;
1506 bestac = newac;
1507 bestpt = newpt;
1508 bestnorm = norm;
1509 }
1510 lastdist = distnew;
1511 lastpt = newpt;
1512 lastac = newac;
1513 move = pt - lastpt;
1514 }
1515 }
1516 }
1517
1518 eval_pt = bestpt;
1519 if (eval_norm)
1520 {
1521 *eval_norm = bestnorm;
1522 }
1523 outside = (nout > 0) ? true : false;
1524 ac = bestac;
1525
1526 return status;
1527 }
1528
1529 //===========================================================================
1530 //Function Name: ac_at_edge
1531 //
1532 //Member Type: PRIVATE
1533 //Description: determine the area coordinate of the facet at the edge
1534 //===========================================================================
ac_at_edge(CartVect & fac,CartVect & eac,int edge_id)1535 void SmoothFace::ac_at_edge(CartVect &fac, // facet area coordinate
1536 CartVect &eac, // edge area coordinate
1537 int edge_id) // id of edge
1538 {
1539 double u, v, w;
1540 switch (edge_id) {
1541 case 0:
1542 u = 0.0;
1543 v = fac[1] / (fac[1] + fac[2]);//v = fac.y() / (fac.y() + fac.z());
1544 w = 1.0 - v;
1545 break;
1546 case 1:
1547 u = fac[0] / (fac[0] + fac[2]);// u = fac.x() / (fac.x() + fac.z());
1548 v = 0.0;
1549 w = 1.0 - u;
1550 break;
1551 case 2:
1552 u = fac[0] / (fac[0] + fac[1]);//u = fac.x() / (fac.x() + fac.y());
1553 v = 1.0 - u;
1554 w = 0.0;
1555 break;
1556 default:
1557 assert(0);
1558 u = -1; // needed to eliminate warnings about used before set
1559 v = -1; // needed to eliminate warnings about used before set
1560 w = -1; // needed to eliminate warnings about used before set
1561 break;
1562 }
1563 eac[0] = u;
1564 eac[1] = v;
1565 eac[2] = w; //= CartVect(u, v, w);
1566 }
1567
1568 //===========================================================================
1569 //Function Name: project_to_facet
1570 //
1571 //Member Type: PUBLIC
1572 //Description: project to a single facet. Uses the input areacoord as
1573 // a starting guess.
1574 //===========================================================================
project_to_facet(EntityHandle facet,CartVect & pt,CartVect & areacoord,CartVect & close_point,bool & outside_facet,double compare_tol)1575 ErrorCode SmoothFace::project_to_facet(EntityHandle facet, CartVect &pt,
1576 CartVect &areacoord, CartVect &close_point, bool &outside_facet,
1577 double compare_tol)
1578 {
1579 const EntityHandle * conn3 = NULL;
1580 int nnodes = 0;
1581 _mb->get_connectivity(facet, conn3, nnodes);
1582 //
1583 //double coords[9]; // store the coordinates for the nodes
1584 //_mb->get_coords(conn3, 3, coords);
1585 CartVect p[3];
1586 _mb->get_coords(conn3, 3, (double*) &p[0]);
1587
1588 int edge_id = -1;
1589 ErrorCode stat = project_to_patch(facet, areacoord, pt, close_point, NULL,
1590 outside_facet, compare_tol, edge_id);
1591 /* }
1592 break;
1593 }*/
1594
1595 return stat;
1596 }
1597
1598 //===========================================================================
1599 //Function Name: is_at_vertex
1600 //
1601 //Member Type: PRIVATE
1602 //Description: determine if the point is at one of the facet's vertices
1603 //===========================================================================
is_at_vertex(EntityHandle facet,CartVect & pt,CartVect & ac,double compare_tol,CartVect & eval_pt,CartVect * eval_norm_ptr)1604 bool SmoothFace::is_at_vertex(EntityHandle facet, // (IN) facet we are evaluating
1605 CartVect &pt, // (IN) the point
1606 CartVect &ac, // (IN) the ac of the point on the facet plane
1607 double compare_tol, // (IN) return TRUE of closer than this
1608 CartVect &eval_pt, // (OUT) location at vertex if TRUE
1609 CartVect *eval_norm_ptr) // (OUT) normal at vertex if TRUE
1610 {
1611 double dist;
1612 CartVect vert_loc;
1613 const double actol = 0.1;
1614 // get coordinates get_coords
1615 const EntityHandle * conn3 = NULL;
1616 int nnodes = 0;
1617 _mb->get_connectivity(facet, conn3, nnodes);
1618 //
1619 //double coords[9]; // store the coordinates for the nodes
1620 //_mb->get_coords(conn3, 3, coords);
1621 CartVect p[3];
1622 _mb->get_coords(conn3, 3, (double*) &p[0]);
1623 // also get the normals at nodes
1624 CartVect NN[3];
1625 _mb->tag_get_data(_gradientTag, conn3, 3, (double*) &NN[0]);
1626 if (fabs(ac[0]) < actol && fabs(ac[1]) < actol)
1627 {
1628 vert_loc = p[2];
1629 dist = (pt - vert_loc).length(); //pt.distance_between( vert_loc );
1630 if (dist <= compare_tol)
1631 {
1632 eval_pt = vert_loc;
1633 if (eval_norm_ptr)
1634 {
1635 *eval_norm_ptr = NN[2];
1636 }
1637 return true;
1638 }
1639 }
1640
1641 if (fabs(ac[0]) < actol && fabs(ac[2]) < actol)
1642 {
1643 vert_loc = p[1];
1644 dist = (pt - vert_loc).length();//pt.distance_between( vert_loc );
1645 if (dist <= compare_tol)
1646 {
1647 eval_pt = vert_loc;
1648 if (eval_norm_ptr)
1649 {
1650 *eval_norm_ptr = NN[1];//facet->point(1)->normal( facet );
1651 }
1652 return true;
1653 }
1654 }
1655
1656 if (fabs(ac[1]) < actol && fabs(ac[2]) < actol)
1657 {
1658 vert_loc = p[0];
1659 dist = (pt - vert_loc).length();//pt.distance_between( vert_loc );
1660 if (dist <= compare_tol)
1661 {
1662 eval_pt = vert_loc;
1663 if (eval_norm_ptr)
1664 {
1665 *eval_norm_ptr = NN[0];
1666 }
1667 return true;
1668 }
1669 }
1670
1671 return false;
1672 }
1673
1674 //===========================================================================
1675 //Function Name: move_ac_inside
1676 //
1677 //Member Type: PRIVATE
1678 //Description: find the closest area coordinate to the boundary of the
1679 // patch if any of its components are < 0
1680 // Return if the ac was modified.
1681 //===========================================================================
move_ac_inside(CartVect & ac,double tol)1682 bool SmoothFace::move_ac_inside(CartVect &ac, double tol)
1683 {
1684 int nout = 0;
1685 if (ac[0] < -tol)
1686 {
1687 ac[0] = 0.0;
1688 ac[1] = ac[1] / (ac[1] + ac[2]); //( ac.y() / (ac.y() + ac.z()) ;
1689 ac[2] = 1. - ac[1]; //ac.z( 1.0 - ac.y() );
1690 nout++;
1691 }
1692 if (ac[1] < -tol)
1693 {
1694 ac[1] = 0.;//ac.y( 0.0 );
1695 ac[0] = ac[0] / (ac[0] + ac[2]);//ac.x( ac.x() / (ac.x() + ac.z()) );
1696 ac[2] = 1. - ac[0];//ac.z( 1.0 - ac.x() );
1697 nout++;
1698 }
1699 if (ac[2] < -tol)
1700 {
1701 ac[2] = 0.;// ac.z( 0.0 );
1702 ac[0] = ac[0] / (ac[0] + ac[1]);//ac.x( ac.x() / (ac.x() + ac.y()) );
1703 ac[1] = 1. - ac[0]; // ac.y( 1.0 - ac.x() );
1704 nout++;
1705 }
1706 return (nout > 0) ? true : false;
1707 }
1708
1709 //===========================================================================
1710 //Function Name: hodograph
1711 //
1712 //Member Type: PUBLIC
1713 //Description: get the hodograph control points for the facet
1714 //Note: This is a triangle cubic patch that is defined by the
1715 // normals of quartic facet control point lattice. Returned coordinates
1716 // in Nijk are defined by the following diagram
1717 //
1718 //
1719 // *9 index polar
1720 // / \ 0 300 point(0)
1721 // / \ 1 210
1722 // 7*-----*8 2 120
1723 // / \ / \ 3 030 point(1)
1724 // / \ / \ 4 201
1725 // 4*----5*-----*6 5 111
1726 // / \ / \ / \ 6 021
1727 // / \ / \ / \ 7 102
1728 // *-----*-----*-----* 8 012
1729 // 0 1 2 3 9 003 point(2)
1730 //
1731
1732 //===========================================================================
1733 //Function Name: eval_bezier_patch_normal
1734 //
1735 //Member Type: PRIVATE
1736 //Description: evaluate the Bezier patch defined at a facet
1737 //===========================================================================
eval_bezier_patch_normal(EntityHandle facet,CartVect & areacoord,CartVect & normal)1738 ErrorCode SmoothFace::eval_bezier_patch_normal(EntityHandle facet,
1739 CartVect &areacoord, CartVect &normal)
1740 {
1741 // interpolate internal control points
1742
1743 CartVect gctrl_pts[6];
1744 //facet->get_control_points( gctrl_pts );
1745 ErrorCode rval = _mb->tag_get_data(_facetCtrlTag, &facet, 1, &gctrl_pts[0]);
1746 assert(rval == MB_SUCCESS);
1747 if (MB_SUCCESS != rval) return rval;
1748 // _gradientTag
1749 // get normals at points
1750 const EntityHandle * conn3 = NULL;
1751 int nnodes = 0;
1752 rval = _mb->get_connectivity(facet, conn3, nnodes);
1753 if (MB_SUCCESS != rval) return rval;
1754
1755 CartVect NN[3];
1756 rval = _mb->tag_get_data(_gradientTag, conn3, 3, &NN[0]);
1757 assert(rval == MB_SUCCESS);
1758 if (MB_SUCCESS != rval) return rval;
1759
1760 if (fabs(areacoord[1] + areacoord[2]) < 1.0e-6)
1761 {
1762 normal = NN[0];
1763 return MB_SUCCESS;
1764 }
1765 if (fabs(areacoord[0] + areacoord[2]) < 1.0e-6)
1766 {
1767 normal = NN[1];//facet->point(1)->normal(facet);
1768 return MB_SUCCESS;
1769 }
1770 if (fabs(areacoord[0] + areacoord[1]) < 1.0e-6)
1771 {
1772 normal = NN[2]; //facet->point(2)->normal(facet);
1773 return MB_SUCCESS;
1774 }
1775
1776 // compute the hodograph of the quartic Gregory patch
1777
1778 CartVect Nijk[10];
1779 //hodograph(facet,areacoord,Nijk);
1780 // start copy from hodograph
1781 //CubitVector gctrl_pts[6];
1782 // facet->get_control_points( gctrl_pts );
1783 CartVect P_facet[3];
1784
1785 //2,1,1
1786 /*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
1787 (areacoord.y() * gctrl_pts[3] +
1788 areacoord.z() * gctrl_pts[4]);*/
1789 P_facet[0] = (1.0e0 / (areacoord[1] + areacoord[2])) * (areacoord[1]
1790 * gctrl_pts[3] + areacoord[2] * gctrl_pts[4]);
1791 //1,2,1
1792 /*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
1793 (areacoord.x() * gctrl_pts[0] +
1794 areacoord.z() * gctrl_pts[5]);*/
1795 P_facet[1] = (1.0e0 / (areacoord[0] + areacoord[2])) * (areacoord[0]
1796 * gctrl_pts[0] + areacoord[2] * gctrl_pts[5]);
1797 //1,1,2
1798 /*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
1799 (areacoord.x() * gctrl_pts[1] +
1800 areacoord.y() * gctrl_pts[2]);*/
1801 P_facet[2] = (1.0e0 / (areacoord[0] + areacoord[1])) * (areacoord[0]
1802 * gctrl_pts[1] + areacoord[1] * gctrl_pts[2]);
1803
1804 // corner control points are just the normals at the points
1805
1806 //3, 0, 0
1807 Nijk[0] = NN[0];
1808 //0, 3, 0
1809 Nijk[3] = NN[1];
1810 //0, 0, 3
1811 Nijk[9] = NN[2];//facet->point(2)->normal(facet);
1812
1813 // fill in the boundary control points. Define as the normal to the local
1814 // triangle formed by the quartic control point lattice
1815
1816 // store here again the 9 control points on the edges
1817 CartVect CP[9]; // 9 control points on the edges,
1818 rval = _mb->tag_get_data(_facetEdgeCtrlTag, &facet, 1, &CP[0]);
1819 if (MB_SUCCESS != rval) return rval;
1820 // there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
1821 //CubitFacetEdge *edge;
1822 //edge = facet->edge( 2 );
1823 //CubitVector ctrl_pts[5];
1824 //edge->control_points(facet, ctrl_pts);
1825
1826 //2, 1, 0
1827 //Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
1828 Nijk[1] = (CP[7] - CP[6]) * (P_facet[0] - CP[6]);
1829 Nijk[1].normalize();
1830
1831 //1, 2, 0
1832 //Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
1833 Nijk[2] = (CP[8] - CP[7]) * (P_facet[1] - CP[7]);
1834 Nijk[2].normalize();
1835
1836 //edge = facet->edge( 0 );
1837 //edge->control_points(facet, ctrl_pts);
1838
1839 //0, 2, 1
1840 //Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
1841 Nijk[6] = (CP[0] - P_facet[1]) * (CP[1] - P_facet[1]);
1842 Nijk[6].normalize();
1843
1844 //0, 1, 2
1845 //Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
1846 Nijk[8] = (CP[1] - P_facet[2]) * (CP[2] - P_facet[2]);
1847 Nijk[8].normalize();
1848
1849 //edge = facet->edge( 1 );
1850 //edge->control_points(facet, ctrl_pts);
1851
1852 //1, 0, 2
1853 //Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
1854 Nijk[7] = (P_facet[2] - CP[4]) * (CP[3] - CP[4]);
1855 Nijk[7].normalize();
1856
1857 //2, 0, 1
1858 //Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
1859 Nijk[4] = (P_facet[0] - CP[5]) * (CP[4] - CP[5]);
1860 Nijk[4].normalize();
1861
1862 //1, 1, 1
1863 Nijk[5] = (P_facet[1] - P_facet[0]) * (P_facet[2] - P_facet[0]);
1864 Nijk[5].normalize();
1865 // end copy from hodograph
1866
1867 // sum the contribution from each of the control points
1868
1869 normal = CartVect(0.0e0, 0.0e0, 0.0e0);
1870
1871 //i=3; j=0; k=0;
1872 //double Bsum = 0.0;
1873 double B = mbcube(areacoord[0]);
1874 //Bsum += B;
1875 normal += B * Nijk[0];
1876
1877 //i=2; j=1; k=0;
1878 B = 3.0 * mbsqr(areacoord[0]) * areacoord[1];
1879 //Bsum += B;
1880 normal += B * Nijk[1];
1881
1882 //i=1; j=2; k=0;
1883 B = 3.0 * areacoord[0] * mbsqr(areacoord[1]);
1884 //Bsum += B;
1885 normal += B * Nijk[2];
1886
1887 //i=0; j=3; k=0;
1888 B = mbcube(areacoord[1]);
1889 //Bsum += B;
1890 normal += B * Nijk[3];
1891
1892 //i=2; j=0; k=1;
1893 B = 3.0 * mbsqr(areacoord[0]) * areacoord[2];
1894 //Bsum += B;
1895 normal += B * Nijk[4];
1896
1897 //i=1; j=1; k=1;
1898 B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
1899 //Bsum += B;
1900 normal += B * Nijk[5];
1901
1902 //i=0; j=2; k=1;
1903 B = 3.0 * mbsqr(areacoord[1]) * areacoord[2];
1904 //Bsum += B;
1905 normal += B * Nijk[6];
1906
1907 //i=1; j=0; k=2;
1908 B = 3.0 * areacoord[0] * mbsqr(areacoord[2]);
1909 //Bsum += B;
1910 normal += B * Nijk[7];
1911
1912 //i=0; j=1; k=2;
1913 B = 3.0 * areacoord[1] * mbsqr(areacoord[2]);
1914 //Bsum += B;
1915 normal += B * Nijk[8];
1916
1917 //i=0; j=0; k=3;
1918 B = mbcube(areacoord[2]);
1919 //Bsum += B;
1920 normal += B * Nijk[9];
1921
1922 //assert(fabs(Bsum - 1.0) < 1e-9);
1923
1924 normal.normalize();
1925
1926 return MB_SUCCESS;
1927 }
1928
get_normals_for_vertices(const EntityHandle * conn2,CartVect N[2])1929 ErrorCode SmoothFace::get_normals_for_vertices(const EntityHandle * conn2,
1930 CartVect N[2])
1931 // this method will be called to retrieve the normals needed in the calculation of control edge points..
1932 {
1933 //CartVect N[2];
1934 //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
1935 ErrorCode rval = _mb->tag_get_data(_gradientTag, conn2, 2, (double*) &N[0]);
1936 return rval;
1937 }
1938
ray_intersection_correct(EntityHandle,CartVect & pt,CartVect & ray,CartVect & eval_pt,double & distance,bool & outside)1939 ErrorCode SmoothFace::ray_intersection_correct(EntityHandle , // (IN) the facet where the patch is defined
1940 CartVect &pt, // (IN) shoot from
1941 CartVect &ray, // (IN) ray direction
1942 CartVect &eval_pt, // (INOUT) The intersection point
1943 double & distance, // (IN OUT) the new distance
1944 bool &outside)
1945 {
1946 // find a point on the smooth surface
1947 CartVect currentPoint = eval_pt;
1948 int numIter = 0;
1949 double improvement = 1.e20;
1950 CartVect diff;
1951 while (numIter++ < 5 && improvement > 0.01)
1952 {
1953 CartVect newPos;
1954
1955 bool trim = false;// is it needed?
1956 outside = true;
1957 CartVect closestPoint;
1958 CartVect normal;
1959
1960 ErrorCode rval = project_to_facets_main(currentPoint, trim, outside,
1961 &newPos, &normal);
1962 if (MB_SUCCESS != rval) return rval;
1963 assert(rval==MB_SUCCESS);
1964 diff = newPos - currentPoint;
1965 improvement = diff.length();
1966 // ( pt + t * ray - closest ) % normal = 0;
1967 // intersect tangent plane that goes through closest point with the direction
1968 // t = normal%(closest-pt) / normal%ray;
1969 double dot = normal % ray; // if it is 0, get out while we can
1970 if (dot < 0.00001)
1971 {
1972 // bad convergence, get out, do not modify anything
1973 return MB_SUCCESS;
1974 }
1975 double t = ((newPos - pt) % normal) / (dot);
1976 currentPoint = pt + t * ray;
1977
1978 }
1979 eval_pt = currentPoint;
1980 diff = currentPoint - pt;
1981 distance = diff.length();
1982 return MB_SUCCESS;
1983 }
1984 }// namespace moab
1985