Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::SmoothCurve Class Reference

#include <SmoothCurve.hpp>

+ Collaboration diagram for moab::SmoothCurve:

Public Member Functions

 SmoothCurve (Interface *mb, EntityHandle curve, GeomTopoTool *gTool)
 
virtual ~SmoothCurve ()
 
virtual double arc_length ()
 
virtual bool is_parametric ()
 Get the parametric status of the curve. More...
 
virtual bool is_periodic (double &period)
 Get the periodic status of the curve. More...
 
virtual void get_param_range (double &u_start, double &u_end)
 Get the parameter range of the curve. More...
 
virtual double u_from_arc_length (double u_root, double arc_length)
 Compute the parameter value at a specified distance along the curve. More...
 
virtual bool position_from_u (double u, double &x, double &y, double &z, double *tg=NULL)
 Evaluate the curve at a specified parameter value. More...
 
virtual void move_to_curve (double &x, double &y, double &z)
 Move a point near the curve to the closest point on the curve. More...
 
virtual double u_from_position (double x, double y, double z, EntityHandle &v, int &indexEdge)
 Get the u parameter value on the curve closest to x,y,z and the point on the curve. More...
 
virtual void start_coordinates (double &x, double &y, double &z)
 Get the starting point of the curve. More...
 
virtual void end_coordinates (double &x, double &y, double &z)
 Get the ending point of the curve. More...
 
void compute_tangents_for_each_edge ()
 
void compute_control_points_on_boundary_edges (double min_dot, std::map< EntityHandle, SmoothFace * > &mapSurfaces, Tag controlPointsTag, Tag markTag)
 
ErrorCode evaluate_smooth_edge (EntityHandle eh, double &tt, CartVect &outv, CartVect &out_tangent)
 

Private Attributes

std::vector< EntityHandle_entities
 
double _leng
 
std::vector< double > _fractions
 
Tag _edgeTag
 
Interface_mb
 
EntityHandle _set
 
GeomTopoTool_gtt
 

Detailed Description

Definition at line 29 of file SmoothCurve.hpp.

Constructor & Destructor Documentation

◆ SmoothCurve()

moab::SmoothCurve::SmoothCurve ( Interface mb,
EntityHandle  curve,
GeomTopoTool gTool 
)

Definition at line 16 of file SmoothCurve.cpp.

17  : _mb( mb ), _set( curve ), _gtt( gTool )
18 {
19  //_mbOut->create_meshset(MESHSET_ORDERED, _oSet);
20  /*_cmlEdgeMesher = new CMLEdgeMesher (this, CML::STANDARD);
21  _cmlEdgeMesher->set_sizing_function(CML::LINEAR_SIZING);*/
22  _leng = 0; // not initialized
23  _edgeTag = 0; // not initialized
24 }

References _edgeTag, and _leng.

◆ ~SmoothCurve()

moab::SmoothCurve::~SmoothCurve ( )
virtual

Definition at line 25 of file SmoothCurve.cpp.

26 {
27  // TODO Auto-generated destructor stub
28 }

Member Function Documentation

◆ arc_length()

double moab::SmoothCurve::arc_length ( )
virtual

Definition at line 30 of file SmoothCurve.cpp.

31 {
32 
33  return _leng;
34 }

References _leng.

◆ compute_control_points_on_boundary_edges()

void moab::SmoothCurve::compute_control_points_on_boundary_edges ( double  min_dot,
std::map< EntityHandle, SmoothFace * > &  mapSurfaces,
Tag  controlPointsTag,
Tag  markTag 
)

Definition at line 428 of file SmoothCurve.cpp.

432 {
433  // these points really need the surfaces they belong to, because the control points on edges
434  // depend on the normals on surfaces
435  // the control points are averaged from different surfaces, by simple mean average
436  // the surfaces have
437  // do we really need min_dot here?
438  // first of all, find out the SmoothFace for each surface set that is adjacent here
439  // GeomTopoTool gTopoTool(_mb);
440  std::vector< EntityHandle > faces;
441  std::vector< int > senses;
442  ErrorCode rval = _gtt->get_senses( _set, faces, senses );
443  if( MB_SUCCESS != rval ) return;
444 
445  // need to find the smooth face attached
446  unsigned int numSurfacesAdjacent = faces.size();
447  // get the edges, and then get the
448  // std::vector<EntityHandle> entities;
449  _mb->get_entities_by_type( _set, MBEDGE, _entities ); // no recursion!!
450  // each edge has the tangent computed already
451  Tag tangentsTag;
452  rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
453  if( rval != MB_SUCCESS ) return; // some error should be thrown
454 
455  // we do not want to search every time
456  std::vector< SmoothFace* > smoothFaceArray;
457  unsigned int i = 0;
458  for( i = 0; i < numSurfacesAdjacent; i++ )
459  {
460  SmoothFace* sms = mapSurfaces[faces[i]];
461  smoothFaceArray.push_back( sms );
462  }
463 
464  unsigned int e = 0;
465  for( e = 0; e < _entities.size(); e++ )
466  {
467  CartVect zero( 0. );
468  CartVect ctrlP[3] = { zero, zero, zero }; // null positions initially
469  // the control points are averaged from connected faces
470  EntityHandle edge = _entities[e]; // the edge in the chain
471 
472  int nnodes;
473  const EntityHandle* conn2;
474  rval = _mb->get_connectivity( edge, conn2, nnodes );
475  if( rval != MB_SUCCESS || 2 != nnodes ) return; // or continue or return error
476 
477  // double coords[6]; // store the coordinates for the nodes
478  CartVect P[2];
479  // ErrorCode rval = _mb->get_coords(conn2, 2, coords);
480  rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
481  if( rval != MB_SUCCESS ) return;
482 
483  CartVect chord = P[1] - P[0];
484  _leng += chord.length();
485  _fractions.push_back( _leng );
486  CartVect N[2];
487 
488  // MBCartVect N0(&normalVec[0]);
489  // MBCartVect N3(&normalVec[3]);
490  CartVect T[2]; // T0, T3
491  // if (edge->num_adj_facets() <= 1) {
492  // stat = compute_curve_tangent(edge, min_dot, T0, T3);
493  // if (stat != CUBIT_SUCCESS)
494  // return stat;
495  //} else {
496  //}
497  rval = _mb->tag_get_data( tangentsTag, &edge, 1, &T[0] );
498  if( rval != MB_SUCCESS ) return;
499 
500  for( i = 0; i < numSurfacesAdjacent; i++ )
501  {
502  CartVect controlForEdge[3];
503  rval = smoothFaceArray[i]->get_normals_for_vertices( conn2, N );
504  if( rval != MB_SUCCESS ) return;
505 
506  rval = smoothFaceArray[i]->init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], controlForEdge );
507  if( rval != MB_SUCCESS ) return;
508 
509  // accumulate those over faces!!!
510  for( int j = 0; j < 3; j++ )
511  {
512  ctrlP[j] += controlForEdge[j];
513  }
514  }
515  // now divide them for the average position!
516  for( int j = 0; j < 3; j++ )
517  {
518  ctrlP[j] /= numSurfacesAdjacent;
519  }
520  // we are done, set the control points now!
521  // edge->control_points(ctrl_pts, 4);
522  rval = _mb->tag_set_data( controlPointsTag, &edge, 1, &ctrlP[0] );
523  if( rval != MB_SUCCESS ) return;
524 
525  this->_edgeTag = controlPointsTag; // this is a tag that will be stored with the edge
526  // is that a waste of memory or not...
527  // also mark the edge for later on
528  unsigned char used = 1;
529  _mb->tag_set_data( markTag, &edge, 1, &used );
530  }
531  // now divide fractions, to make them vary from 0 to 1
532  assert( _leng > 0. );
533  for( e = 0; e < _entities.size(); e++ )
534  _fractions[e] /= _leng;
535 }

References _edgeTag, _entities, _fractions, _gtt, _leng, _mb, _set, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), moab::GeomTopoTool::get_senses(), moab::CartVect::length(), MB_SUCCESS, MB_TYPE_DOUBLE, MBEDGE, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::Interface::tag_set_data().

Referenced by moab::FBEngine::initializeSmoothing().

◆ compute_tangents_for_each_edge()

void moab::SmoothCurve::compute_tangents_for_each_edge ( )

Definition at line 377 of file SmoothCurve.cpp.

378 {
379  // will retrieve the edges in each set; they are retrieved in order they were put into, because
380  // these sets are "MESHSET_ORDERED"
381  // retrieve the tag handle for the tangents; it should have been created already
382  // this tangents are computed for the chain of edges that form a geometric edge
383  // some checks should be performed on the vertices, but we trust the correctness of the model
384  // completely (like the vertices should match in the chain...)
385  Tag tangentsTag;
386  ErrorCode rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
387  if( rval != MB_SUCCESS ) return; // some error should be thrown
388  std::vector< EntityHandle > entities;
389  _mb->get_entities_by_type( _set, MBEDGE, entities ); // no recursion!!
390  // basically, each tangent at a node will depend on previous tangent
391  int nbEdges = entities.size();
392  // now, we can advance in the loop
393  // the only special problem is if the first node coincides with the last node, then we should
394  // consider the closed loop; or maybe we should look at angles in that case too?
395  // also, should we look at the 2 semi-circles case? How to decide if we need to continue the
396  // "tangents" maybe we can do that later, and we can alter the tangents at the feature nodes, in
397  // the directions of the loops again, do we need to decide the "closed" loop or not? Not yet...
398  EntityHandle previousEdge = entities[0]; // this is the first edge in the chain
399  CartVect TP[2]; // tangents for the previous edge
400  rval = _mb->tag_get_data( tangentsTag, &previousEdge, 1, &TP[0] ); // tangents for previous edge
401  if( rval != MB_SUCCESS ) return; // some error should be thrown
402  CartVect TC[2]; // tangents for the current edge
403  EntityHandle currentEdge;
404  for( int i = 1; i < nbEdges; i++ )
405  {
406  // current edge will start after first one
407  currentEdge = entities[i];
408  rval = _mb->tag_get_data( tangentsTag, &currentEdge, 1, &TC[0] ); //
409  if( rval != MB_SUCCESS ) return; // some error should be thrown
410  // now compute the new tangent at common vertex; reset tangents for previous edge and
411  // current edge a little bit of CPU and memory waste, but this is life
412  CartVect T = 0.5 * TC[0] + 0.5 * TP[1]; //
413  T.normalize();
414  TP[1] = T;
415  rval = _mb->tag_set_data( tangentsTag, &previousEdge, 1, &TP[0] ); //
416  if( rval != MB_SUCCESS ) return; // some error should be thrown
417  TC[0] = T;
418  rval = _mb->tag_set_data( tangentsTag, &currentEdge, 1, &TC[0] ); //
419  if( rval != MB_SUCCESS ) return; // some error should be thrown
420  // now set the next edge
421  previousEdge = currentEdge;
422  TP[0] = TC[0];
423  TP[1] = TC[1];
424  }
425  return;
426 }

References _mb, _set, entities, ErrorCode, moab::Interface::get_entities_by_type(), MB_SUCCESS, MB_TYPE_DOUBLE, MBEDGE, moab::CartVect::normalize(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and TC.

Referenced by moab::FBEngine::initializeSmoothing().

◆ end_coordinates()

void moab::SmoothCurve::end_coordinates ( double &  x,
double &  y,
double &  z 
)
virtual

Get the ending point of the curve.

Parameters
xThe x coordinate of the start point
yThe y coordinate of the start point
zThe z coordinate of the start point

Definition at line 360 of file SmoothCurve.cpp.

361 {
362 
363  int nnodes = 0;
364  const EntityHandle* conn2 = NULL;
365  _mb->get_connectivity( _entities[_entities.size() - 1], conn2, nnodes );
366  double c[3];
367  // careful, the second node here
368  _mb->get_coords( &conn2[1], 1, c );
369 
370  x = c[0];
371  y = c[1];
372  z = c[2];
373  return;
374 }

References _entities, _mb, moab::Interface::get_connectivity(), and moab::Interface::get_coords().

◆ evaluate_smooth_edge()

ErrorCode moab::SmoothCurve::evaluate_smooth_edge ( EntityHandle  eh,
double &  tt,
CartVect outv,
CartVect out_tangent 
)

Definition at line 537 of file SmoothCurve.cpp.

538 {
539  CartVect P[2]; // P0 and P1
540  CartVect controlPoints[3]; // edge control points
541  double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
542 
543  // project the position to the linear edge
544  // t is from 0 to 1 only!!
545  // double tt = (t + 1) * 0.5;
546  if( tt <= 0.0 ) tt = 0.0;
547  if( tt >= 1.0 ) tt = 1.0;
548 
549  int nnodes = 0;
550  const EntityHandle* conn2 = NULL;
551  ErrorCode rval = _mb->get_connectivity( eh, conn2, nnodes );
552  if( rval != MB_SUCCESS ) return rval;
553 
554  rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
555  if( rval != MB_SUCCESS ) return rval;
556 
557  if( 0 == _edgeTag )
558  {
559  rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
560  if( rval != MB_SUCCESS ) return rval;
561  }
562  rval = _mb->tag_get_data( _edgeTag, &eh, 1, (double*)&controlPoints[0] );
563  if( rval != MB_SUCCESS ) return rval;
564 
565  t2 = tt * tt;
566  t3 = t2 * tt;
567  t4 = t3 * tt;
568  one_minus_t = 1. - tt;
569  one_minus_t2 = one_minus_t * one_minus_t;
570  one_minus_t3 = one_minus_t2 * one_minus_t;
571  one_minus_t4 = one_minus_t3 * one_minus_t;
572 
573  outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
574  4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
575 
576  out_tangent = -4. * one_minus_t3 * P[0] + 4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
577  12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
578  4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
579  return MB_SUCCESS;
580 }

References _edgeTag, _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), MB_SUCCESS, MB_TYPE_DOUBLE, moab::Interface::tag_get_data(), and moab::Interface::tag_get_handle().

Referenced by position_from_u().

◆ get_param_range()

void moab::SmoothCurve::get_param_range ( double &  u_start,
double &  u_end 
)
virtual

Get the parameter range of the curve.

Parameters
u_startThe beginning curve parameter
u_endThe ending curve parameter
Note
The numerical value of u_start may be greater than the numerical value of u_end.

Definition at line 70 of file SmoothCurve.cpp.

71 {
72  // assert(_ref_edge);
73  u_start = 0;
74  u_end = 1.;
75 
76  return;
77 }

Referenced by moab::FBEngine::getEntURange().

◆ is_parametric()

bool moab::SmoothCurve::is_parametric ( )
virtual

Get the parametric status of the curve.

Returns
true if curve is parametric, false otherwise.

Definition at line 39 of file SmoothCurve.cpp.

40 {
41  return true;
42 }

◆ is_periodic()

bool moab::SmoothCurve::is_periodic ( double &  period)
virtual

Get the periodic status of the curve.

Parameters
periodThe period of the curve if periodic.
Returns
true if curve is periodic, false otherwise.

Definition at line 49 of file SmoothCurve.cpp.

50 {
51  // assert(_ref_edge);
52  // return _ref_edge->is_periodic( period);
53  Range vsets;
54  _mb->get_child_meshsets( _set, vsets ); // num_hops =1
55  if( vsets.size() == 1 )
56  {
57  period = _leng;
58  return true; // true , especially for ice sheet data
59  }
60  return false;
61 }

References _leng, _mb, _set, moab::Interface::get_child_meshsets(), and moab::Range::size().

◆ move_to_curve()

void moab::SmoothCurve::move_to_curve ( double &  x,
double &  y,
double &  z 
)
virtual

Move a point near the curve to the closest point on the curve.

Parameters
xThe x coordinate of the point
yThe y coordinate of the point
zThe z coordinate of the point

Definition at line 144 of file SmoothCurve.cpp.

145 {
146 
147  // find closest point to the curve, and the parametric position
148  // must be close by, but how close ???
149  EntityHandle v;
150  int edgeIndex;
151  double u = u_from_position( x, y, z, v, edgeIndex );
152  position_from_u( u, x, y, z );
153 
154  return;
155 }

References position_from_u(), and u_from_position().

Referenced by moab::FBEngine::getEntClosestPt().

◆ position_from_u()

bool moab::SmoothCurve::position_from_u ( double  u,
double &  x,
double &  y,
double &  z,
double *  tg = NULL 
)
virtual

Evaluate the curve at a specified parameter value.

Parameters
uThe parameter at which to evaluate the curve
xThe x coordinate of the evaluated point
yThe y coordinate of the evaluated point
zThe z coordinate of the evaluated point

Definition at line 105 of file SmoothCurve.cpp.

106 {
107 
108  // _fractions are increasing, so find the
109  double* ptr = std::lower_bound( &_fractions[0], ( &_fractions[0] ) + _fractions.size(), u );
110  int index = ptr - &_fractions[0];
111  double nextFraction = _fractions[index];
112  double prevFraction = 0;
113  if( index > 0 )
114  {
115  prevFraction = _fractions[index - 1];
116  }
117  double t = ( u - prevFraction ) / ( nextFraction - prevFraction );
118 
119  EntityHandle edge = _entities[index];
120 
121  CartVect position, tangent;
122  ErrorCode rval = evaluate_smooth_edge( edge, t, position, tangent );
123  if( MB_SUCCESS != rval ) return false;
124  assert( rval == MB_SUCCESS );
125  x = position[0];
126  y = position[1];
127  z = position[2];
128  if( tg )
129  {
130  // we need to do some scaling,
131  double dtdu = 1 / ( nextFraction - prevFraction );
132  tg[0] = tangent[0] * dtdu;
133  tg[1] = tangent[1] * dtdu;
134  tg[2] = tangent[2] * dtdu;
135  }
136 
137  return true;
138 }

References _entities, _fractions, ErrorCode, evaluate_smooth_edge(), and MB_SUCCESS.

Referenced by moab::FBEngine::getEntTgntU(), moab::FBEngine::getEntUtoXYZ(), and move_to_curve().

◆ start_coordinates()

void moab::SmoothCurve::start_coordinates ( double &  x,
double &  y,
double &  z 
)
virtual

Get the starting point of the curve.

Parameters
xThe x coordinate of the start point
yThe y coordinate of the start point
zThe z coordinate of the start point

Definition at line 339 of file SmoothCurve.cpp.

340 {
341 
342  int nnodes = 0;
343  const EntityHandle* conn2 = NULL;
344  _mb->get_connectivity( _entities[0], conn2, nnodes );
345  double c[3];
346  _mb->get_coords( conn2, 1, c );
347 
348  x = c[0];
349  y = c[1];
350  z = c[2];
351 
352  return;
353 }

References _entities, _mb, moab::Interface::get_connectivity(), and moab::Interface::get_coords().

◆ u_from_arc_length()

double moab::SmoothCurve::u_from_arc_length ( double  u_root,
double  arc_leng 
)
virtual

Compute the parameter value at a specified distance along the curve.

Parameters
u_rootThe start parameter from which to compute the distance along the curve.
arc_lengthThe distance to move along the curve.
Note
For positive values of arc_length the distance will be computed in the direction of increasing parameter value along the curve. For negative values of arc_length the distance will be computed in the direction of decreasing parameter value along the curve.
Returns
The parametric coordinate u along the curve

Definition at line 92 of file SmoothCurve.cpp.

93 {
94 
95  if( _leng <= 0 ) return 0;
96  return u_root + arc_leng / _leng;
97 }

References _leng.

◆ u_from_position()

double moab::SmoothCurve::u_from_position ( double  x,
double  y,
double  z,
EntityHandle v,
int &  edgeIndex 
)
virtual

Get the u parameter value on the curve closest to x,y,z and the point on the curve.

Parameters
xThe x coordinate of the point
yThe y coordinate of the point
zThe z coordinate of the point
Returns
The parametric coordinate u on the curve

Definition at line 165 of file SmoothCurve.cpp.

166 {
167  // this is an iterative process, expensive usually
168  // get first all nodes , and their positions
169  // find the closest node (and edge), and from there do some
170  // iterations up to a point
171  // do not exaggerate with convergence criteria
172 
173  v = 0; // we do not have a close by vertex yet
174  CartVect initialPos( x, y, z );
175  double u = 0;
176  int nbNodes = (int)_entities.size() * 2; // the mesh edges are stored
177  std::vector< EntityHandle > nodesConnec;
178  nodesConnec.resize( nbNodes );
179  ErrorCode rval = this->_mb->get_connectivity( &( _entities[0] ), nbNodes / 2, nodesConnec );
180  if( MB_SUCCESS != rval )
181  {
182  std::cout << "error in getting connectivity\n";
183  return 0;
184  }
185  // collapse nodesConnec, nodes should be in order
186  for( int k = 0; k < nbNodes / 2; k++ )
187  {
188  nodesConnec[k + 1] = nodesConnec[2 * k + 1];
189  }
190  int numNodes = nbNodes / 2 + 1;
191  std::vector< CartVect > coordNodes;
192  coordNodes.resize( numNodes );
193 
194  rval = _mb->get_coords( &( nodesConnec[0] ), numNodes, (double*)&( coordNodes[0] ) );
195  if( MB_SUCCESS != rval )
196  {
197  std::cout << "error in getting node positions\n";
198  return 0;
199  }
200  // find the closest node, then find the closest edge, based on closest node
201 
202  int indexNode = 0;
203  double minDist = 1.e30;
204  // expensive linear search
205  for( int i = 0; i < numNodes; i++ )
206  {
207  double d1 = ( initialPos - coordNodes[i] ).length();
208  if( d1 < minDist )
209  {
210  indexNode = i;
211  minDist = d1;
212  }
213  }
214  double tolerance = 0.00001; // what is the unit?
215  // something reasonable
216  if( minDist < tolerance )
217  {
218  v = nodesConnec[indexNode];
219  // we are done, just return the proper u (from fractions)
220  if( indexNode == 0 )
221  {
222  return 0; // first node has u = 0
223  }
224  else
225  return _fractions[indexNode - 1]; // fractions[0] > 0!!)
226  }
227  // find the mesh edge; could be previous or next edge
228  edgeIndex = indexNode; // could be the previous one, though!!
229  if( edgeIndex == numNodes - 1 )
230  edgeIndex--; // we have one less edge, and do not worry about next edge
231  else
232  {
233  if( edgeIndex > 0 )
234  {
235  // could be the previous; decide based on distance to the other
236  // nodes of the 2 connected edges
237  CartVect prevNodePos = coordNodes[edgeIndex - 1];
238  CartVect nextNodePos = coordNodes[edgeIndex + 1];
239  if( ( prevNodePos - initialPos ).length_squared() < ( nextNodePos - initialPos ).length_squared() )
240  {
241  edgeIndex--;
242  }
243  }
244  }
245  // now, we know for sure that the closest point is somewhere on edgeIndex edge
246  //
247 
248  // do newton iteration for local t between 0 and 1
249 
250  // copy from evaluation method
251  CartVect P[2]; // P0 and P1
252  CartVect controlPoints[3]; // edge control points
253  double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
254 
255  P[0] = coordNodes[edgeIndex];
256  P[1] = coordNodes[edgeIndex + 1];
257 
258  if( 0 == _edgeTag )
259  {
260  rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
261  if( rval != MB_SUCCESS ) return 0;
262  }
263  rval = _mb->tag_get_data( _edgeTag, &( _entities[edgeIndex] ), 1, (double*)&controlPoints[0] );
264  if( rval != MB_SUCCESS ) return rval;
265 
266  // starting point
267  double tt = 0.5; // between the 2 ends of the edge
268  int iterations = 0;
269  // find iteratively a better point
270  int maxIterations = 10; // not too many
271  CartVect outv;
272  // we will solve minimize F = 0.5 * ( ini - r(t) )^2
273  // so solve F'(t) = 0
274  // Iteration: t_ -> t - F'(t)/F"(t)
275  // F'(t) = r'(t) (ini-r(t) )
276  // F"(t) = r"(t) (ini-r(t) ) - (r'(t))^2
277  while( iterations < maxIterations )
278  //
279  {
280  t2 = tt * tt;
281  t3 = t2 * tt;
282  t4 = t3 * tt;
283  one_minus_t = 1. - tt;
284  one_minus_t2 = one_minus_t * one_minus_t;
285  one_minus_t3 = one_minus_t2 * one_minus_t;
286  one_minus_t4 = one_minus_t3 * one_minus_t;
287 
288  outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] +
289  6. * one_minus_t2 * t2 * controlPoints[1] + 4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
290 
291  CartVect out_tangent = -4. * one_minus_t3 * P[0] +
292  4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
293  12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
294  4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
295 
296  CartVect second_deriv =
297  12. * one_minus_t2 * P[0] +
298  4. * ( -3. * one_minus_t2 - 3. * one_minus_t2 + 6. * tt * one_minus_t ) * controlPoints[0] +
299  12. * ( one_minus_t2 - 4 * tt * one_minus_t + t2 ) * controlPoints[1] +
300  4. * ( 6. * tt - 12 * t2 ) * controlPoints[2] + 12. * t2 * P[1];
301  CartVect diff = outv - initialPos;
302  double F_d = out_tangent % diff;
303  double F_dd = second_deriv % diff + out_tangent.length_squared();
304 
305  if( 0 == F_dd ) break; // get out, we found minimum?
306 
307  double delta_t = -F_d / F_dd;
308 
309  if( fabs( delta_t ) < 0.000001 ) break;
310  tt = tt + delta_t;
311  if( tt < 0 )
312  {
313  tt = 0.;
314  v = nodesConnec[edgeIndex]; // we are at end of mesh edge
315  break;
316  }
317  if( tt > 1 )
318  {
319  tt = 1;
320  v = nodesConnec[edgeIndex + 1]; // we are at one end
321  break;
322  }
323  iterations++;
324  }
325  // so we have t on the segment, convert to u, which should
326  // be between _fractions[edgeIndex] numbers
327  double prevFraction = 0;
328  if( edgeIndex > 0 ) prevFraction = _fractions[edgeIndex - 1];
329 
330  u = prevFraction + tt * ( _fractions[edgeIndex] - prevFraction );
331  return u;
332 }

References _edgeTag, _entities, _fractions, _mb, ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), length(), moab::CartVect::length_squared(), length_squared(), MB_SUCCESS, MB_TYPE_DOUBLE, moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), and moab::tolerance.

Referenced by move_to_curve(), and moab::FBEngine::split_edge_at_point().

Member Data Documentation

◆ _edgeTag

Tag moab::SmoothCurve::_edgeTag
private

◆ _entities

std::vector< EntityHandle > moab::SmoothCurve::_entities
private

◆ _fractions

std::vector< double > moab::SmoothCurve::_fractions
private

◆ _gtt

GeomTopoTool* moab::SmoothCurve::_gtt
private

Definition at line 138 of file SmoothCurve.hpp.

Referenced by compute_control_points_on_boundary_edges().

◆ _leng

double moab::SmoothCurve::_leng
private

◆ _mb

◆ _set

EntityHandle moab::SmoothCurve::_set
private

The documentation for this class was generated from the following files: