Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SmoothCurve.cpp
Go to the documentation of this file.
1 /*
2  * SmoothCurve.cpp
3  *
4  */
5 
6 #include "SmoothCurve.hpp"
7 //#include "SmoothVertex.hpp"
8 #include "SmoothFace.hpp"
9 #include <cassert>
10 #include <iostream>
11 #include "moab/GeomTopoTool.hpp"
12 
13 namespace moab
14 {
15 
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 }
26 {
27  // TODO Auto-generated destructor stub
28 }
29 
31 {
32 
33  return _leng;
34 }
35 
36 //! \brief Get the parametric status of the curve.
37 //!
38 //! \return \a true if curve is parametric, \a false otherwise.
40 {
41  return true;
42 }
43 
44 //! \brief Get the periodic status of the curve.
45 //!
46 //! \param period The period of the curve if periodic.
47 //!
48 //! \return \a true if curve is periodic, \a false otherwise.
49 bool SmoothCurve::is_periodic( double& period )
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 }
62 
63 //! \brief Get the parameter range of the curve.
64 //!
65 //! \param u_start The beginning curve parameter
66 //! \param u_end The ending curve parameter
67 //!
68 //! \note The numerical value of \a u_start may be greater
69 //! than the numerical value of \a u_end.
70 void SmoothCurve::get_param_range( double& u_start, double& u_end )
71 {
72  // assert(_ref_edge);
73  u_start = 0;
74  u_end = 1.;
75 
76  return;
77 }
78 
79 //! Compute the parameter value at a specified distance along the curve.
80 //!
81 //! \param u_root The start parameter from which to compute the distance
82 //! along the curve.
83 //! \param arc_length The distance to move along the curve.
84 //!
85 //! \note For positive values of \a arc_length the distance will be
86 //! computed in the direction of increasing parameter value along the
87 //! curve. For negative values of \a arc_length the distance will be
88 //! computed in the direction of decreasing parameter value along the
89 //! curve.
90 //!
91 //! \return The parametric coordinate u along the curve
92 double SmoothCurve::u_from_arc_length( double u_root, double arc_leng )
93 {
94 
95  if( _leng <= 0 ) return 0;
96  return u_root + arc_leng / _leng;
97 }
98 
99 //! \brief Evaluate the curve at a specified parameter value.
100 //!
101 //! \param u The parameter at which to evaluate the curve
102 //! \param x The x coordinate of the evaluated point
103 //! \param y The y coordinate of the evaluated point
104 //! \param z The z coordinate of the evaluated point
105 bool SmoothCurve::position_from_u( double u, double& x, double& y, double& z, double* tg )
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 }
139 //! \brief Move a point near the curve to the closest point on the curve.
140 //!
141 //! \param x The x coordinate of the point
142 //! \param y The y coordinate of the point
143 //! \param z The z coordinate of the point
144 void SmoothCurve::move_to_curve( double& x, double& y, double& z )
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 }
156 
157 //! Get the u parameter value on the curve closest to x,y,z
158 //! and the point on the curve.
159 //!
160 //! \param x The x coordinate of the point
161 //! \param y The y coordinate of the point
162 //! \param z The z coordinate of the point
163 //!
164 //! \return The parametric coordinate u on the curve
165 double SmoothCurve::u_from_position( double x, double y, double z, EntityHandle& v, int& edgeIndex )
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 }
333 
334 //! \brief Get the starting point of the curve.
335 //!
336 //! \param x The x coordinate of the start point
337 //! \param y The y coordinate of the start point
338 //! \param z The z coordinate of the start point
339 void SmoothCurve::start_coordinates( double& x, double& y, double& z )
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 }
354 
355 //! \brief Get the ending point of the curve.
356 //!
357 //! \param x The x coordinate of the start point
358 //! \param y The y coordinate of the start point
359 //! \param z The z coordinate of the start point
360 void SmoothCurve::end_coordinates( double& x, double& y, double& z )
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 }
375 
376 // this will recompute the 2 tangents for each edge, considering the geo edge they are into
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 }
427 
429  std::map< EntityHandle, SmoothFace* >& mapSurfaces,
430  Tag controlPointsTag,
431  Tag markTag )
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 }
536 
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 }
581 } // namespace moab