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