Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SmoothFace.cpp
Go to the documentation of this file.
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