Branch data Line data Source code
1 : : /*
2 : : * SmoothCurve.cpp
3 : : *
4 : : */
5 : :
6 : : #include "SmoothCurve.hpp"
7 : : //#include "SmoothVertex.hpp"
8 : : #include "SmoothFace.hpp"
9 : : #include "assert.h"
10 : : #include <iostream>
11 : : #include "moab/GeomTopoTool.hpp"
12 : :
13 : : namespace moab
14 : : {
15 : :
16 : 30 : SmoothCurve::SmoothCurve( Interface* mb, EntityHandle curve, GeomTopoTool* gTool )
17 [ + - ]: 30 : : _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 : 30 : _leng = 0; // not initialized
23 : 30 : _edgeTag = 0; // not initialized
24 : 30 : }
25 : 90 : SmoothCurve::~SmoothCurve()
26 : : {
27 : : // TODO Auto-generated destructor stub
28 [ - + ]: 60 : }
29 : :
30 : 0 : double SmoothCurve::arc_length()
31 : : {
32 : :
33 : 0 : 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 : 0 : bool SmoothCurve::is_parametric()
40 : : {
41 : 0 : 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 : 0 : bool SmoothCurve::is_periodic( double& period )
50 : : {
51 : : // assert(_ref_edge);
52 : : // return _ref_edge->is_periodic( period);
53 [ # # ]: 0 : Range vsets;
54 [ # # ]: 0 : _mb->get_child_meshsets( _set, vsets ); // num_hops =1
55 [ # # ][ # # ]: 0 : if( vsets.size() == 1 )
56 : : {
57 : 0 : period = _leng;
58 : 0 : return true; // true , especially for ice sheet data
59 : : }
60 : 0 : 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 : 0 : void SmoothCurve::get_param_range( double& u_start, double& u_end )
71 : : {
72 : : // assert(_ref_edge);
73 : 0 : u_start = 0;
74 : 0 : u_end = 1.;
75 : :
76 : 0 : 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 : 0 : double SmoothCurve::u_from_arc_length( double u_root, double arc_leng )
93 : : {
94 : :
95 [ # # ]: 0 : if( _leng <= 0 ) return 0;
96 : 0 : 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 : 24 : 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 [ + - ][ + - ]: 24 : double* ptr = std::lower_bound( &_fractions[0], ( &_fractions[0] ) + _fractions.size(), u );
[ + - ]
110 [ + - ]: 24 : int index = ptr - &_fractions[0];
111 [ + - ]: 24 : double nextFraction = _fractions[index];
112 : 24 : double prevFraction = 0;
113 [ + - ][ + - ]: 24 : if( index > 0 ) { prevFraction = _fractions[index - 1]; }
114 : 24 : double t = ( u - prevFraction ) / ( nextFraction - prevFraction );
115 : :
116 [ + - ]: 24 : EntityHandle edge = _entities[index];
117 : :
118 [ + - ][ + - ]: 24 : CartVect position, tangent;
119 [ + - ]: 24 : ErrorCode rval = evaluate_smooth_edge( edge, t, position, tangent );
120 [ - + ]: 24 : if( MB_SUCCESS != rval ) return false;
121 [ - + ]: 24 : assert( rval == MB_SUCCESS );
122 [ + - ]: 24 : x = position[0];
123 [ + - ]: 24 : y = position[1];
124 [ + - ]: 24 : z = position[2];
125 [ - + ]: 24 : if( tg )
126 : : {
127 : : // we need to do some scaling,
128 : 0 : double dtdu = 1 / ( nextFraction - prevFraction );
129 [ # # ]: 0 : tg[0] = tangent[0] * dtdu;
130 [ # # ]: 0 : tg[1] = tangent[1] * dtdu;
131 [ # # ]: 0 : tg[2] = tangent[2] * dtdu;
132 : : }
133 : :
134 : 24 : return true;
135 : : }
136 : : //! \brief Move a point near the curve to the closest point on the curve.
137 : : //!
138 : : //! \param x The x coordinate of the point
139 : : //! \param y The y coordinate of the point
140 : : //! \param z The z coordinate of the point
141 : 12 : void SmoothCurve::move_to_curve( double& x, double& y, double& z )
142 : : {
143 : :
144 : : // find closest point to the curve, and the parametric position
145 : : // must be close by, but how close ???
146 : : EntityHandle v;
147 : : int edgeIndex;
148 [ + - ]: 12 : double u = u_from_position( x, y, z, v, edgeIndex );
149 [ + - ]: 12 : position_from_u( u, x, y, z );
150 : :
151 : 12 : return;
152 : : }
153 : :
154 : : //! Get the u parameter value on the curve closest to x,y,z
155 : : //! and the point on the curve.
156 : : //!
157 : : //! \param x The x coordinate of the point
158 : : //! \param y The y coordinate of the point
159 : : //! \param z The z coordinate of the point
160 : : //!
161 : : //! \return The parametric coordinate u on the curve
162 : 12 : double SmoothCurve::u_from_position( double x, double y, double z, EntityHandle& v, int& edgeIndex )
163 : : {
164 : : // this is an iterative process, expensive usually
165 : : // get first all nodes , and their positions
166 : : // find the closest node (and edge), and from there do some
167 : : // iterations up to a point
168 : : // do not exaggerate with convergence criteria
169 : :
170 : 12 : v = 0; // we do not have a close by vertex yet
171 [ + - ]: 12 : CartVect initialPos( x, y, z );
172 : 12 : double u = 0;
173 : 12 : int nbNodes = (int)_entities.size() * 2; // the mesh edges are stored
174 [ + - ]: 12 : std::vector< EntityHandle > nodesConnec;
175 [ + - ]: 12 : nodesConnec.resize( nbNodes );
176 [ + - ][ + - ]: 12 : ErrorCode rval = this->_mb->get_connectivity( &( _entities[0] ), nbNodes / 2, nodesConnec );
177 [ - + ]: 12 : if( MB_SUCCESS != rval )
178 : : {
179 [ # # ]: 0 : std::cout << "error in getting connectivity\n";
180 : 0 : return 0;
181 : : }
182 : : // collapse nodesConnec, nodes should be in order
183 [ + + ]: 530 : for( int k = 0; k < nbNodes / 2; k++ )
184 : : {
185 [ + - ][ + - ]: 518 : nodesConnec[k + 1] = nodesConnec[2 * k + 1];
186 : : }
187 : 12 : int numNodes = nbNodes / 2 + 1;
188 [ + - ]: 24 : std::vector< CartVect > coordNodes;
189 [ + - ]: 12 : coordNodes.resize( numNodes );
190 : :
191 [ + - ][ + - ]: 12 : rval = _mb->get_coords( &( nodesConnec[0] ), numNodes, (double*)&( coordNodes[0] ) );
[ + - ]
192 [ - + ]: 12 : if( MB_SUCCESS != rval )
193 : : {
194 [ # # ]: 0 : std::cout << "error in getting node positions\n";
195 : 0 : return 0;
196 : : }
197 : : // find the closest node, then find the closest edge, based on closest node
198 : :
199 : 12 : int indexNode = 0;
200 : 12 : double minDist = 1.e30;
201 : : // expensive linear search
202 [ + + ]: 542 : for( int i = 0; i < numNodes; i++ )
203 : : {
204 [ + - ][ + - ]: 530 : double d1 = ( initialPos - coordNodes[i] ).length();
[ + - ]
205 [ + + ]: 530 : if( d1 < minDist )
206 : : {
207 : 176 : indexNode = i;
208 : 176 : minDist = d1;
209 : : }
210 : : }
211 : 12 : double tolerance = 0.00001; // what is the unit?
212 : : // something reasonable
213 [ - + ]: 12 : if( minDist < tolerance )
214 : : {
215 [ # # ]: 0 : v = nodesConnec[indexNode];
216 : : // we are done, just return the proper u (from fractions)
217 [ # # ]: 0 : if( indexNode == 0 )
218 : : {
219 : 0 : return 0; // first node has u = 0
220 : : }
221 : : else
222 [ # # ]: 0 : return _fractions[indexNode - 1]; // fractions[0] > 0!!)
223 : : }
224 : : // find the mesh edge; could be previous or next edge
225 : 12 : edgeIndex = indexNode; // could be the previous one, though!!
226 [ - + ]: 12 : if( edgeIndex == numNodes - 1 )
227 : 0 : edgeIndex--; // we have one less edge, and do not worry about next edge
228 : : else
229 : : {
230 [ + - ]: 12 : if( edgeIndex > 0 )
231 : : {
232 : : // could be the previous; decide based on distance to the other
233 : : // nodes of the 2 connected edges
234 [ + - ]: 12 : CartVect prevNodePos = coordNodes[edgeIndex - 1];
235 [ + - ]: 12 : CartVect nextNodePos = coordNodes[edgeIndex + 1];
236 [ + - ][ + - ]: 12 : if( ( prevNodePos - initialPos ).length_squared() < ( nextNodePos - initialPos ).length_squared() )
[ + - ][ + - ]
[ + + ]
237 : 12 : { edgeIndex--; }
238 : : }
239 : : }
240 : : // now, we know for sure that the closest point is somewhere on edgeIndex edge
241 : : //
242 : :
243 : : // do newton iteration for local t between 0 and 1
244 : :
245 : : // copy from evaluation method
246 [ + - ][ + + ]: 36 : CartVect P[2]; // P0 and P1
247 [ + - ][ + + ]: 48 : CartVect controlPoints[3]; // edge control points
248 : : double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
249 : :
250 [ + - ]: 12 : P[0] = coordNodes[edgeIndex];
251 [ + - ]: 12 : P[1] = coordNodes[edgeIndex + 1];
252 : :
253 [ - + ]: 12 : if( 0 == _edgeTag )
254 : : {
255 [ # # ]: 0 : rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
256 [ # # ]: 0 : if( rval != MB_SUCCESS ) return 0;
257 : : }
258 [ + - ][ + - ]: 12 : rval = _mb->tag_get_data( _edgeTag, &( _entities[edgeIndex] ), 1, (double*)&controlPoints[0] );
259 [ - + ]: 12 : if( rval != MB_SUCCESS ) return rval;
260 : :
261 : : // starting point
262 : 12 : double tt = 0.5; // between the 2 ends of the edge
263 : 12 : int iterations = 0;
264 : : // find iteratively a better point
265 : 12 : int maxIterations = 10; // not too many
266 [ + - ]: 12 : CartVect outv;
267 : : // we will solve minimize F = 0.5 * ( ini - r(t) )^2
268 : : // so solve F'(t) = 0
269 : : // Iteration: t_ -> t - F'(t)/F"(t)
270 : : // F'(t) = r'(t) (ini-r(t) )
271 : : // F"(t) = r"(t) (ini-r(t) ) - (r'(t))^2
272 [ + - ]: 31 : while( iterations < maxIterations )
273 : : //
274 : : {
275 : 31 : t2 = tt * tt;
276 : 31 : t3 = t2 * tt;
277 : 31 : t4 = t3 * tt;
278 : 31 : one_minus_t = 1. - tt;
279 : 31 : one_minus_t2 = one_minus_t * one_minus_t;
280 : 31 : one_minus_t3 = one_minus_t2 * one_minus_t;
281 : 31 : one_minus_t4 = one_minus_t3 * one_minus_t;
282 : :
283 [ + - ][ + - ]: 31 : outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] +
[ + - ][ + - ]
284 [ + - ][ + - ]: 62 : 6. * one_minus_t2 * t2 * controlPoints[1] + 4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
[ + - ][ + - ]
[ + - ]
285 : :
286 [ + - ][ + - ]: 31 : CartVect out_tangent = -4. * one_minus_t3 * P[0] +
287 [ + - ][ + - ]: 62 : 4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
288 [ + - ][ + - ]: 62 : 12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
289 [ + - ][ + - ]: 62 : 4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
[ + - ]
290 : :
291 : : CartVect second_deriv =
292 [ + - ][ + - ]: 31 : 12. * one_minus_t2 * P[0] +
293 [ + - ][ + - ]: 62 : 4. * ( -3. * one_minus_t2 - 3. * one_minus_t2 + 6. * tt * one_minus_t ) * controlPoints[0] +
294 [ + - ][ + - ]: 62 : 12. * ( one_minus_t2 - 4 * tt * one_minus_t + t2 ) * controlPoints[1] +
295 [ + - ][ + - ]: 62 : 4. * ( 6. * tt - 12 * t2 ) * controlPoints[2] + 12. * t2 * P[1];
[ + - ]
296 [ + - ]: 31 : CartVect diff = outv - initialPos;
297 [ + - ]: 31 : double F_d = out_tangent % diff;
298 [ + - ][ + - ]: 31 : double F_dd = second_deriv % diff + out_tangent.length_squared();
299 : :
300 [ - + ]: 43 : if( 0 == F_dd ) break; // get out, we found minimum?
301 : :
302 : 31 : double delta_t = -F_d / F_dd;
303 : :
304 [ + + ]: 31 : if( fabs( delta_t ) < 0.000001 ) break;
305 : 21 : tt = tt + delta_t;
306 [ + + ]: 21 : if( tt < 0 )
307 : : {
308 : 2 : tt = 0.;
309 [ + - ]: 2 : v = nodesConnec[edgeIndex]; // we are at end of mesh edge
310 : 2 : break;
311 : : }
312 [ - + ]: 19 : if( tt > 1 )
313 : : {
314 : 0 : tt = 1;
315 [ # # ]: 0 : v = nodesConnec[edgeIndex + 1]; // we are at one end
316 : 0 : break;
317 : : }
318 : 19 : iterations++;
319 : : }
320 : : // so we have t on the segment, convert to u, which should
321 : : // be between _fractions[edgeIndex] numbers
322 : 12 : double prevFraction = 0;
323 [ + - ][ + - ]: 12 : if( edgeIndex > 0 ) prevFraction = _fractions[edgeIndex - 1];
324 : :
325 [ + - ]: 12 : u = prevFraction + tt * ( _fractions[edgeIndex] - prevFraction );
326 : 24 : return u;
327 : : }
328 : :
329 : : //! \brief Get the starting point of the curve.
330 : : //!
331 : : //! \param x The x coordinate of the start point
332 : : //! \param y The y coordinate of the start point
333 : : //! \param z The z coordinate of the start point
334 : 0 : void SmoothCurve::start_coordinates( double& x, double& y, double& z )
335 : : {
336 : :
337 : 0 : int nnodes = 0;
338 : 0 : const EntityHandle* conn2 = NULL;
339 [ # # ][ # # ]: 0 : _mb->get_connectivity( _entities[0], conn2, nnodes );
340 : : double c[3];
341 [ # # ]: 0 : _mb->get_coords( conn2, 1, c );
342 : :
343 : 0 : x = c[0];
344 : 0 : y = c[1];
345 : 0 : z = c[2];
346 : :
347 : 0 : return;
348 : : }
349 : :
350 : : //! \brief Get the ending point of the curve.
351 : : //!
352 : : //! \param x The x coordinate of the start point
353 : : //! \param y The y coordinate of the start point
354 : : //! \param z The z coordinate of the start point
355 : 0 : void SmoothCurve::end_coordinates( double& x, double& y, double& z )
356 : : {
357 : :
358 : 0 : int nnodes = 0;
359 : 0 : const EntityHandle* conn2 = NULL;
360 [ # # ][ # # ]: 0 : _mb->get_connectivity( _entities[_entities.size() - 1], conn2, nnodes );
361 : : double c[3];
362 : : // careful, the second node here
363 [ # # ]: 0 : _mb->get_coords( &conn2[1], 1, c );
364 : :
365 : 0 : x = c[0];
366 : 0 : y = c[1];
367 : 0 : z = c[2];
368 : 0 : return;
369 : : }
370 : :
371 : : // this will recompute the 2 tangents for each edge, considering the geo edge they are into
372 : 30 : void SmoothCurve::compute_tangents_for_each_edge()
373 : : {
374 : : // will retrieve the edges in each set; they are retrieved in order they were put into, because
375 : : // these sets are "MESHSET_ORDERED"
376 : : // retrieve the tag handle for the tangents; it should have been created already
377 : : // this tangents are computed for the chain of edges that form a geometric edge
378 : : // some checks should be performed on the vertices, but we trust the correctness of the model
379 : : // completely (like the vertices should match in the chain...)
380 : : Tag tangentsTag;
381 [ + - ]: 30 : ErrorCode rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
382 [ - + ]: 30 : if( rval != MB_SUCCESS ) return; // some error should be thrown
383 [ + - ]: 30 : std::vector< EntityHandle > entities;
384 [ + - ]: 30 : _mb->get_entities_by_type( _set, MBEDGE, entities ); // no recursion!!
385 : : // basically, each tangent at a node will depend on previous tangent
386 : 30 : int nbEdges = entities.size();
387 : : // now, we can advance in the loop
388 : : // the only special problem is if the first node coincides with the last node, then we should
389 : : // consider the closed loop; or maybe we should look at angles in that case too?
390 : : // also, should we look at the 2 semi-circles case? How to decide if we need to continue the
391 : : // "tangents" maybe we can do that later, and we can alter the tangents at the feature nodes, in
392 : : // the directions of the loops again, do we need to decide the "closed" loop or not? Not yet...
393 [ + - ]: 30 : EntityHandle previousEdge = entities[0]; // this is the first edge in the chain
394 [ + - ][ + + ]: 90 : CartVect TP[2]; // tangents for the previous edge
395 [ + - ]: 30 : rval = _mb->tag_get_data( tangentsTag, &previousEdge, 1, &TP[0] ); // tangents for previous edge
396 [ - + ]: 30 : if( rval != MB_SUCCESS ) return; // some error should be thrown
397 [ + - ][ + + ]: 90 : CartVect TC[2]; // tangents for the current edge
398 : : EntityHandle currentEdge;
399 [ + + ]: 2651 : for( int i = 1; i < nbEdges; i++ )
400 : : {
401 : : // current edge will start after first one
402 [ + - ]: 2621 : currentEdge = entities[i];
403 [ + - ]: 2621 : rval = _mb->tag_get_data( tangentsTag, ¤tEdge, 1, &TC[0] ); //
404 [ - + ]: 2621 : if( rval != MB_SUCCESS ) return; // some error should be thrown
405 : : // now compute the new tangent at common vertex; reset tangents for previous edge and
406 : : // current edge a little bit of CPU and memory waste, but this is life
407 [ + - ][ + - ]: 2621 : CartVect T = 0.5 * TC[0] + 0.5 * TP[1]; //
[ + - ]
408 [ + - ]: 2621 : T.normalize();
409 : 2621 : TP[1] = T;
410 [ + - ]: 2621 : rval = _mb->tag_set_data( tangentsTag, &previousEdge, 1, &TP[0] ); //
411 [ - + ]: 2621 : if( rval != MB_SUCCESS ) return; // some error should be thrown
412 : 2621 : TC[0] = T;
413 [ + - ]: 2621 : rval = _mb->tag_set_data( tangentsTag, ¤tEdge, 1, &TC[0] ); //
414 [ - + ]: 2621 : if( rval != MB_SUCCESS ) return; // some error should be thrown
415 : : // now set the next edge
416 : 2621 : previousEdge = currentEdge;
417 : 2621 : TP[0] = TC[0];
418 : 2621 : TP[1] = TC[1];
419 : : }
420 : 30 : return;
421 : : }
422 : :
423 : 30 : void SmoothCurve::compute_control_points_on_boundary_edges( double, std::map< EntityHandle, SmoothFace* >& mapSurfaces,
424 : : Tag controlPointsTag, Tag markTag )
425 : : {
426 : : // these points really need the surfaces they belong to, because the control points on edges
427 : : // depend on the normals on surfaces
428 : : // the control points are averaged from different surfaces, by simple mean average
429 : : // the surfaces have
430 : : // do we really need min_dot here?
431 : : // first of all, find out the SmoothFace for each surface set that is adjacent here
432 : : // GeomTopoTool gTopoTool(_mb);
433 [ + - ]: 30 : std::vector< EntityHandle > faces;
434 [ + - ][ + - ]: 60 : std::vector< int > senses;
435 [ + - ]: 30 : ErrorCode rval = _gtt->get_senses( _set, faces, senses );
436 [ - + ]: 30 : if( MB_SUCCESS != rval ) return;
437 : :
438 : : // need to find the smooth face attached
439 : 30 : unsigned int numSurfacesAdjacent = faces.size();
440 : : // get the edges, and then get the
441 : : // std::vector<EntityHandle> entities;
442 [ + - ]: 30 : _mb->get_entities_by_type( _set, MBEDGE, _entities ); // no recursion!!
443 : : // each edge has the tangent computed already
444 : : Tag tangentsTag;
445 [ + - ]: 30 : rval = _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, tangentsTag );
446 [ - + ]: 30 : if( rval != MB_SUCCESS ) return; // some error should be thrown
447 : :
448 : : // we do not want to search every time
449 [ + - ][ + - ]: 60 : std::vector< SmoothFace* > smoothFaceArray;
450 : 30 : unsigned int i = 0;
451 [ + + ]: 63 : for( i = 0; i < numSurfacesAdjacent; i++ )
452 : : {
453 [ + - ][ + - ]: 33 : SmoothFace* sms = mapSurfaces[faces[i]];
454 [ + - ]: 33 : smoothFaceArray.push_back( sms );
455 : : }
456 : :
457 : 30 : unsigned int e = 0;
458 [ + + ]: 2681 : for( e = 0; e < _entities.size(); e++ )
459 : : {
460 [ + - ]: 2651 : CartVect zero( 0. );
461 : 2651 : CartVect ctrlP[3] = { zero, zero, zero }; // null positions initially
462 : : // the control points are averaged from connected faces
463 [ + - ]: 2651 : EntityHandle edge = _entities[e]; // the edge in the chain
464 : :
465 : : int nnodes;
466 : : const EntityHandle* conn2;
467 [ + - ]: 2651 : rval = _mb->get_connectivity( edge, conn2, nnodes );
468 [ + - ][ - + ]: 2651 : if( rval != MB_SUCCESS || 2 != nnodes ) return; // or continue or return error
469 : :
470 : : // double coords[6]; // store the coordinates for the nodes
471 [ + - ][ + + ]: 7953 : CartVect P[2];
472 : : // ErrorCode rval = _mb->get_coords(conn2, 2, coords);
473 [ + - ]: 2651 : rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
474 [ - + ]: 2651 : if( rval != MB_SUCCESS ) return;
475 : :
476 [ + - ]: 2651 : CartVect chord = P[1] - P[0];
477 [ + - ]: 2651 : _leng += chord.length();
478 [ + - ]: 2651 : _fractions.push_back( _leng );
479 [ + - ][ + + ]: 7953 : CartVect N[2];
480 : :
481 : : // MBCartVect N0(&normalVec[0]);
482 : : // MBCartVect N3(&normalVec[3]);
483 [ + - ][ + + ]: 7953 : CartVect T[2]; // T0, T3
484 : : // if (edge->num_adj_facets() <= 1) {
485 : : // stat = compute_curve_tangent(edge, min_dot, T0, T3);
486 : : // if (stat != CUBIT_SUCCESS)
487 : : // return stat;
488 : : //} else {
489 : : //}
490 [ + - ]: 2651 : rval = _mb->tag_get_data( tangentsTag, &edge, 1, &T[0] );
491 [ - + ]: 2651 : if( rval != MB_SUCCESS ) return;
492 : :
493 [ + + ]: 5318 : for( i = 0; i < numSurfacesAdjacent; i++ )
494 : : {
495 [ + - ][ + + ]: 10668 : CartVect controlForEdge[3];
496 [ + - ][ + - ]: 2667 : rval = smoothFaceArray[i]->get_normals_for_vertices( conn2, N );
497 [ - + ]: 2667 : if( rval != MB_SUCCESS ) return;
498 : :
499 [ + - ][ + - ]: 2667 : rval = smoothFaceArray[i]->init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], controlForEdge );
500 [ - + ]: 2667 : if( rval != MB_SUCCESS ) return;
501 : :
502 : : // accumulate those over faces!!!
503 [ + + ]: 10668 : for( int j = 0; j < 3; j++ )
504 : : {
505 [ + - ]: 8001 : ctrlP[j] += controlForEdge[j];
506 : : }
507 : : }
508 : : // now divide them for the average position!
509 [ + + ]: 10604 : for( int j = 0; j < 3; j++ )
510 : : {
511 [ + - ]: 7953 : ctrlP[j] /= numSurfacesAdjacent;
512 : : }
513 : : // we are done, set the control points now!
514 : : // edge->control_points(ctrl_pts, 4);
515 [ + - ]: 2651 : rval = _mb->tag_set_data( controlPointsTag, &edge, 1, &ctrlP[0] );
516 [ - + ]: 2651 : if( rval != MB_SUCCESS ) return;
517 : :
518 : 2651 : this->_edgeTag = controlPointsTag; // this is a tag that will be stored with the edge
519 : : // is that a waste of memory or not...
520 : : // also mark the edge for later on
521 : 2651 : unsigned char used = 1;
522 [ + - ]: 2651 : _mb->tag_set_data( markTag, &edge, 1, &used );
523 : : }
524 : : // now divide fractions, to make them vary from 0 to 1
525 [ - + ]: 30 : assert( _leng > 0. );
526 [ + + ][ + - ]: 2681 : for( e = 0; e < _entities.size(); e++ )
527 [ + - ]: 2681 : _fractions[e] /= _leng;
528 : : }
529 : :
530 : 24 : ErrorCode SmoothCurve::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv, CartVect& out_tangent )
531 : : {
532 [ + - ][ + + ]: 72 : CartVect P[2]; // P0 and P1
533 [ + - ][ + + ]: 96 : CartVect controlPoints[3]; // edge control points
534 : : double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
535 : :
536 : : // project the position to the linear edge
537 : : // t is from 0 to 1 only!!
538 : : // double tt = (t + 1) * 0.5;
539 [ - + ]: 24 : if( tt <= 0.0 ) tt = 0.0;
540 [ + + ]: 24 : if( tt >= 1.0 ) tt = 1.0;
541 : :
542 : 24 : int nnodes = 0;
543 : 24 : const EntityHandle* conn2 = NULL;
544 [ + - ]: 24 : ErrorCode rval = _mb->get_connectivity( eh, conn2, nnodes );
545 [ - + ]: 24 : if( rval != MB_SUCCESS ) return rval;
546 : :
547 [ + - ]: 24 : rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
548 [ - + ]: 24 : if( rval != MB_SUCCESS ) return rval;
549 : :
550 [ - + ]: 24 : if( 0 == _edgeTag )
551 : : {
552 [ # # ]: 0 : rval = _mb->tag_get_handle( "CONTROLEDGE", 9, MB_TYPE_DOUBLE, _edgeTag );
553 [ # # ]: 0 : if( rval != MB_SUCCESS ) return rval;
554 : : }
555 [ + - ]: 24 : rval = _mb->tag_get_data( _edgeTag, &eh, 1, (double*)&controlPoints[0] );
556 [ - + ]: 24 : if( rval != MB_SUCCESS ) return rval;
557 : :
558 : 24 : t2 = tt * tt;
559 : 24 : t3 = t2 * tt;
560 : 24 : t4 = t3 * tt;
561 : 24 : one_minus_t = 1. - tt;
562 : 24 : one_minus_t2 = one_minus_t * one_minus_t;
563 : 24 : one_minus_t3 = one_minus_t2 * one_minus_t;
564 : 24 : one_minus_t4 = one_minus_t3 * one_minus_t;
565 : :
566 [ + - ][ + - ]: 24 : outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
[ + - ][ + - ]
[ + - ][ + - ]
567 [ + - ][ + - ]: 48 : 4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
[ + - ]
568 : :
569 [ + - ][ + - ]: 24 : out_tangent = -4. * one_minus_t3 * P[0] + 4. * ( one_minus_t3 - 3. * tt * one_minus_t2 ) * controlPoints[0] +
[ + - ][ + - ]
570 [ + - ][ + - ]: 48 : 12. * ( tt * one_minus_t2 - t2 * one_minus_t ) * controlPoints[1] +
571 [ + - ][ + - ]: 48 : 4. * ( 3. * t2 * one_minus_t - t3 ) * controlPoints[2] + 4. * t3 * P[1];
[ + - ]
572 : 24 : return MB_SUCCESS;
573 : : }
574 [ + - ][ + - ]: 44 : } // namespace moab
|
