 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*
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
16 SmoothCurve::SmoothCurve( Interface* mb, EntityHandle curve, GeomTopoTool* gTool )
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 }
25 SmoothCurve::~SmoothCurve()
26 {
27 // TODO Auto-generated destructor stub
28 }
29
30 double SmoothCurve::arc_length()
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.
39 bool SmoothCurve::is_parametric()
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
377 void SmoothCurve::compute_tangents_for_each_edge()
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, ¤tEdge, 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, ¤tEdge, 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
428 void SmoothCurve::compute_control_points_on_boundary_edges( double,
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
537 ErrorCode SmoothCurve::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv, CartVect& out_tangent )
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
1.9.1.