Branch data Line data Source code
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 : 8040 : bool within_tolerance( CartVect& p1, CartVect& p2, const double& tolerance )
19 : : {
20 [ + + ]: 8040 : if( ( fabs( p1[0] - p2[0] ) < tolerance ) && ( fabs( p1[1] - p2[1] ) < tolerance ) &&
[ - + # # ]
[ - + ]
21 : 0 : ( fabs( p1[2] - p2[2] ) < tolerance ) )
22 : 0 : return true;
23 : 8040 : return false;
24 : : }
25 : 0 : int numAdjTriInSet( Interface* mb, EntityHandle startEdge, EntityHandle set )
26 : : {
27 [ # # ]: 0 : std::vector< EntityHandle > adjTri;
28 [ # # ]: 0 : mb->get_adjacencies( &startEdge, 1, 2, false, adjTri, Interface::UNION );
29 : : // count how many are in the set
30 : 0 : int nbInSet = 0;
31 [ # # ]: 0 : for( size_t i = 0; i < adjTri.size(); i++ )
32 : : {
33 [ # # ]: 0 : EntityHandle tri = adjTri[i];
34 [ # # ][ # # ]: 0 : if( mb->contains_entities( set, &tri, 1 ) ) nbInSet++;
35 : : }
36 : 0 : return nbInSet;
37 : : }
38 : :
39 : : bool debug_surf_eval1 = false;
40 : :
41 : 12 : 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 [ + - ][ + - ]: 12 : _evaluationsCounter( 0 )
45 : : {
46 : : //_smooth_face = NULL;
47 : : //_mbOut->create_meshset(MESHSET_SET, _oSet); //will contain the
48 : : // get also the obb_root
49 [ + - ]: 12 : if( _my_geomTopoTool )
50 : : {
51 [ + - ]: 12 : _my_geomTopoTool->get_root( this->_set, _obb_root );
52 [ - + ][ # # ]: 12 : if( debug_surf_eval1 ) _my_geomTopoTool->obb_tree()->stats( _obb_root, std::cout );
[ # # ]
53 : : }
54 : 12 : }
55 : :
56 [ - + ]: 48 : SmoothFace::~SmoothFace() {}
57 : :
58 : 0 : double SmoothFace::area()
59 : : {
60 : : // find the area of this entity
61 : : // assert(_smooth_face);
62 : : // double area1 = _smooth_face->area();
63 : 0 : double totArea = 0.;
64 [ # # ][ # # ]: 0 : for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
[ # # ][ # # ]
[ # # ]
65 : : {
66 [ # # ]: 0 : EntityHandle tria = *it;
67 : : const EntityHandle* conn3;
68 : : int nnodes;
69 [ # # ]: 0 : _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 [ # # ][ # # ]: 0 : CartVect p[3];
74 [ # # ]: 0 : _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 [ # # ]: 0 : CartVect AB( p[1] - p[0] ); //(p2 - p1);
80 [ # # ]: 0 : CartVect BC( p[2] - p[1] ); //(p3 - p2);
81 [ # # ]: 0 : CartVect normal = AB * BC;
82 [ # # ]: 0 : totArea += normal.length() / 2;
83 : : }
84 : 0 : return totArea;
85 : : }
86 : : // these tags will be collected for deletion
87 : 6 : void SmoothFace::append_smooth_tags( std::vector< Tag >& smoothTags )
88 : : {
89 : : // these are created locally, for each smooth face
90 : 6 : smoothTags.push_back( _gradientTag );
91 : 6 : smoothTags.push_back( _planeTag );
92 : 6 : }
93 : 0 : void SmoothFace::bounding_box( double box_min[3], double box_max[3] )
94 : : {
95 : :
96 [ # # ]: 0 : for( int i = 0; i < 3; i++ )
97 : : {
98 : 0 : box_min[i] = _minim[i];
99 : 0 : box_max[i] = _maxim[i];
100 : : }
101 : : // _minim, _maxim
102 : 0 : }
103 : :
104 : 1763 : void SmoothFace::move_to_surface( double& x, double& y, double& z )
105 : : {
106 : :
107 [ + - ]: 1763 : CartVect loc2( x, y, z );
108 : 1763 : bool trim = false; // is it needed?
109 : 1763 : bool outside = true;
110 [ + - ]: 1763 : CartVect closestPoint;
111 : :
112 [ + - ]: 1763 : ErrorCode rval = project_to_facets_main( loc2, trim, outside, &closestPoint, NULL );
113 [ - + ]: 3526 : if( MB_SUCCESS != rval ) return;
114 [ - + ]: 1763 : assert( rval == MB_SUCCESS );
115 [ + - ]: 1763 : x = closestPoint[0];
116 [ + - ]: 1763 : y = closestPoint[1];
117 [ + - ]: 1763 : 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 : 8 : bool SmoothFace::normal_at( double x, double y, double z, double& nx, double& ny, double& nz )
128 : : {
129 : :
130 [ + - ]: 8 : CartVect loc2( x, y, z );
131 : :
132 : 8 : bool trim = false; // is it needed?
133 : 8 : bool outside = true;
134 : : // CartVect closestPoint;// not needed
135 : : // not interested in normal
136 [ + - ]: 8 : CartVect normal;
137 [ + - ]: 8 : ErrorCode rval = project_to_facets_main( loc2, trim, outside, NULL, &normal );
138 [ - + ]: 8 : if( MB_SUCCESS != rval ) return false;
139 [ - + ]: 8 : assert( rval == MB_SUCCESS );
140 [ + - ]: 8 : nx = normal[0];
141 [ + - ]: 8 : ny = normal[1];
142 [ + - ]: 8 : nz = normal[2];
143 : :
144 : 8 : return true;
145 : : }
146 : :
147 : 12 : ErrorCode SmoothFace::compute_control_points_on_edges( double min_dot, Tag edgeCtrlTag, Tag markTag )
148 : : {
149 : :
150 : 12 : _edgeCtrlTag = edgeCtrlTag;
151 : 12 : _markTag = markTag;
152 : :
153 : : // now, compute control points for all edges that are not marked already (they are no on the
154 : : // boundary!)
155 [ + - ][ + - ]: 68205 : for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
[ + - ][ + - ]
[ + + ]
156 : : {
157 [ + - ]: 68193 : EntityHandle edg = *it;
158 : : // is the edge marked? already computed
159 : 68193 : unsigned char tagVal = 0;
160 [ + - ]: 68193 : _mb->tag_get_data( _markTag, &edg, 1, &tagVal );
161 [ + + ]: 68193 : if( tagVal ) continue;
162 : : // double min_dot;
163 [ + - ]: 65526 : init_bezier_edge( edg, min_dot );
164 : 65526 : tagVal = 1;
165 [ + - ]: 65526 : _mb->tag_set_data( _markTag, &edg, 1, &tagVal );
166 : : }
167 : 12 : return MB_SUCCESS;
168 : : }
169 : :
170 : 12 : 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 [ - + ]: 12 : if( NULL == _mb ) return 1; // fail
176 [ + - ]: 12 : _triangles.clear();
177 [ + - ]: 12 : ErrorCode rval = _mb->get_entities_by_type( _set, MBTRI, _triangles );
178 [ - + ]: 12 : if( MB_SUCCESS != rval ) return 1;
179 : : // get a new range of edges, and decide the loops from here
180 [ + - ]: 12 : _edges.clear();
181 [ + - ]: 12 : rval = _mb->get_adjacencies( _triangles, 1, true, _edges, Interface::UNION );
182 [ - + ]: 12 : assert( rval == MB_SUCCESS );
183 : :
184 [ + - ]: 12 : rval = _mb->get_adjacencies( _triangles, 0, false, _nodes, Interface::UNION );
185 [ - + ]: 12 : assert( rval == MB_SUCCESS );
186 : :
187 : : // initialize bounding box limits
188 [ + - ]: 12 : CartVect vert1;
189 [ + - ]: 12 : EntityHandle v1 = _nodes[0]; // first vertex
190 [ + - ]: 12 : _mb->get_coords( &v1, 1, (double*)&vert1 );
191 : 12 : _minim = vert1;
192 : 12 : _maxim = vert1;
193 : :
194 : 12 : 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 [ + - ]: 12 : unsigned long setId = _mb->id_from_handle( _set );
199 : 12 : char name[50] = { 0 };
200 : : sprintf( name, "GRADIENT%lu",
201 : 12 : setId ); // name should be something like GRADIENT29, where 29 is the set ID of the face
202 [ + - ]: 12 : rval = _mb->tag_get_handle( name, 3, MB_TYPE_DOUBLE, _gradientTag, MB_TAG_DENSE | MB_TAG_CREAT, &defNormal );
203 [ - + ]: 12 : assert( rval == MB_SUCCESS );
204 : :
205 : 12 : double defPlane[4] = { 0., 0., 1., 0. };
206 : : // also define a plane tag ; this will be for each triangle
207 : 12 : char namePlaneTag[50] = { 0 };
208 : 12 : sprintf( namePlaneTag, "PLANE%lu", setId );
209 [ + - ]: 12 : rval = _mb->tag_get_handle( "PLANE", 4, MB_TYPE_DOUBLE, _planeTag, MB_TAG_DENSE | MB_TAG_CREAT, &defPlane );
210 [ - + ]: 12 : 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 [ + - ][ + - ]: 44585 : for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
[ + - ][ + - ]
[ + + ]
214 : : {
215 [ + - ]: 44573 : EntityHandle tria = *it;
216 : : const EntityHandle* conn3;
217 : : int nnodes;
218 [ + - ]: 44573 : _mb->get_connectivity( tria, conn3, nnodes );
219 [ - + ]: 44573 : if( nnodes != 3 ) return 1; // error
220 : : // double coords[9]; // store the coordinates for the nodes
221 : : //_mb->get_coords(conn3, 3, coords);
222 [ + - ][ + + ]: 178292 : CartVect p[3];
223 [ + - ]: 44573 : _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 [ + - ]: 44573 : CartVect AB( p[1] - p[0] ); //(p2 - p1);
229 [ + - ]: 44573 : CartVect BC( p[2] - p[1] ); //(p3 - p2);
230 [ + - ]: 44573 : CartVect CA( p[0] - p[2] ); //(p1 - p3);
231 : : double a[3];
232 [ + - ][ + - ]: 44573 : a[1] = angle( AB, -BC ); // angle at B (p2), etc.
233 [ + - ][ + - ]: 44573 : a[2] = angle( BC, -CA );
234 [ + - ][ + - ]: 44573 : a[0] = angle( CA, -AB );
235 [ + - ][ + - ]: 44573 : CartVect normal = -AB * CA;
236 [ + - ]: 44573 : normal.normalize();
237 : : double plane[4];
238 : :
239 [ + - ]: 44573 : const double* coordNormal = normal.array();
240 : :
241 : 44573 : plane[0] = coordNormal[0];
242 : 44573 : plane[1] = coordNormal[1];
243 : 44573 : plane[2] = coordNormal[2];
244 [ + - ][ + - ]: 44573 : plane[3] = -normal % p[0]; // dot product
245 : : // set the plane
246 [ + - ]: 44573 : rval = _mb->tag_set_data( _planeTag, &tria, 1, plane );
247 [ - + ]: 44573 : 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 [ + - ]: 44573 : _mb->tag_get_data( _gradientTag, conn3, 3, values );
253 [ + + ]: 178292 : 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 : 133719 : values[3 * i + 0] += a[i] * coordNormal[0];
257 : 133719 : values[3 * i + 1] += a[i] * coordNormal[1];
258 : 133719 : values[3 * i + 2] += a[i] * coordNormal[2];
259 : : }
260 : :
261 : : // reset those values
262 [ + - ]: 44573 : _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 [ + - ]: 12 : int numNodes = _nodes.size();
267 [ + - ][ + - ]: 12 : double* normalVal = new double[3 * numNodes];
268 : : _mb->tag_get_data( _gradientTag, _nodes,
269 [ + - ]: 12 : normalVal ); // get all the normal values at the _nodes
270 [ + + ]: 23644 : for( int i = 0; i < numNodes; i++ )
271 : : {
272 [ + - ]: 23632 : CartVect p1( &normalVal[3 * i] );
273 [ + - ]: 23632 : p1.normalize();
274 [ + - ]: 23632 : p1.get( &normalVal[3 * i] );
275 : : }
276 : :
277 : : // reset the normal values after normalization
278 [ + - ]: 12 : _mb->tag_set_data( _gradientTag, _nodes, normalVal );
279 : : // print the loops size and some other stuff
280 [ - + ]: 12 : if( debug_surf_eval1 )
281 : : {
282 [ # # ][ # # ]: 0 : std::cout << " normals at " << numNodes << " nodes" << std::endl;
[ # # ][ # # ]
283 : 0 : int i = 0;
284 [ # # ][ # # ]: 0 : for( Range::iterator it = _nodes.begin(); it != _nodes.end(); ++it, i++ )
[ # # ][ # # ]
[ # # ]
285 : : {
286 [ # # ]: 0 : EntityHandle node = *it;
287 [ # # ][ # # ]: 0 : std::cout << " Node id " << _mb->id_from_handle( node ) << " " << normalVal[3 * i] << " "
[ # # ][ # # ]
[ # # ][ # # ]
288 [ # # ][ # # ]: 0 : << normalVal[3 * i + 1] << " " << normalVal[3 * i + 2] << std::endl;
[ # # ][ # # ]
289 : : }
290 : : }
291 : :
292 [ + - ]: 12 : delete[] normalVal;
293 : :
294 : 12 : 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 : 65526 : 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 [ + - ][ + + ]: 262104 : CartVect ctrl_pts[3];
313 : 65526 : int nnodes = 0;
314 : 65526 : const EntityHandle* conn2 = NULL;
315 [ + - ]: 65526 : ErrorCode rval = _mb->get_connectivity( edge, conn2, nnodes );
316 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
317 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
318 : :
319 [ - + ]: 65526 : assert( 2 == nnodes );
320 : : // double coords[6]; // store the coordinates for the nodes
321 [ + - ][ + + ]: 196578 : CartVect P[2];
322 : : // ErrorCode rval = _mb->get_coords(conn2, 2, coords);
323 [ + - ]: 65526 : rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
324 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
325 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
326 : :
327 : : // CartVect P0(&coords[0]);
328 : : // CartVect P3(&coords[3]);
329 : :
330 : : // double normalVec[6];
331 [ + - ][ + + ]: 196578 : CartVect N[2];
332 : : //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
333 [ + - ]: 65526 : rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
334 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
335 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
336 : :
337 [ + - ][ + + ]: 196578 : CartVect T[2]; // T0, T3
338 : :
339 [ + - ]: 65526 : rval = _mb->tag_get_data( _tangentsTag, &edge, 1, &T[0] );
340 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
341 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
342 : :
343 [ + - ]: 65526 : rval = init_edge_control_points( P[0], P[1], N[0], N[1], T[0], T[1], ctrl_pts );
344 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
345 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
346 : :
347 [ + - ]: 65526 : rval = _mb->tag_set_data( _edgeCtrlTag, &edge, 1, &ctrl_pts[0] );
348 [ - + ]: 65526 : assert( rval == MB_SUCCESS );
349 [ - + ]: 65526 : if( MB_SUCCESS != rval ) return rval;
350 : :
351 [ - + ]: 65526 : if( debug_surf_eval1 )
352 : : {
353 [ # # ][ # # ]: 0 : std::cout << "edge: " << _mb->id_from_handle( edge ) << " tangents: " << T[0] << T[1] << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
354 [ # # ][ # # ]: 0 : std::cout << " points: " << P[0] << " " << P[1] << std::endl;
[ # # ][ # # ]
[ # # ]
355 [ # # ][ # # ]: 0 : std::cout << " normals: " << N[0] << " " << N[1] << std::endl;
[ # # ][ # # ]
[ # # ]
356 [ # # ][ # # ]: 0 : std::cout << " Control points " << ctrl_pts[0] << " " << ctrl_pts[1] << " " << ctrl_pts[2] << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
357 : : }
358 : :
359 : 65526 : return MB_SUCCESS;
360 : : }
361 : :
362 : 12 : ErrorCode SmoothFace::compute_tangents_for_each_edge()
363 : : // they will be used for control points
364 : : {
365 : 12 : double defTangents[6] = { 0., 0., 0., 0., 0., 0. };
366 : : ErrorCode rval =
367 [ + - ]: 12 : _mb->tag_get_handle( "TANGENTS", 6, MB_TYPE_DOUBLE, _tangentsTag, MB_TAG_DENSE | MB_TAG_CREAT, &defTangents );
368 [ - + ]: 12 : 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 [ + - ][ + - ]: 68205 : for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
[ + - ][ + - ]
[ + + ]
372 : : {
373 [ + - ]: 68193 : EntityHandle edg = *it;
374 : :
375 : : int nnodes;
376 : : const EntityHandle* conn2; //
377 [ + - ]: 68193 : _mb->get_connectivity( edg, conn2, nnodes );
378 [ - + ]: 68193 : assert( nnodes == 2 );
379 [ + - ][ + + ]: 204579 : CartVect P[2]; // store the coordinates for the nodes
380 [ + - ]: 68193 : rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
381 [ - + ]: 68193 : if( MB_SUCCESS != rval ) return rval;
382 [ - + ]: 68193 : assert( rval == MB_SUCCESS );
383 [ + - ][ + + ]: 204579 : CartVect T[2];
384 [ + - ]: 68193 : T[0] = P[1] - P[0];
385 [ + - ]: 68193 : T[0].normalize();
386 : 68193 : T[1] = T[0]; //
387 : : _mb->tag_set_data( _tangentsTag, &edg, 1,
388 [ + - ]: 68193 : (double*)&T[0] ); // set the tangents computed at every edge
389 : : }
390 : 12 : 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 : 68193 : ErrorCode SmoothFace::init_edge_control_points( CartVect& P0, CartVect& P3, CartVect& N0, CartVect& N3, CartVect& T0,
400 : : CartVect& T3, CartVect* Pi )
401 : : {
402 [ + - ][ + + ]: 340965 : CartVect Vi[4];
403 : 68193 : Vi[0] = P0;
404 : 68193 : Vi[3] = P3;
405 [ + - ]: 68193 : CartVect P03( P3 - P0 );
406 [ + - ]: 68193 : double di = P03.length();
407 [ + - ]: 68193 : double ai = N0 % N3; // this is the dot operator, the same as in cgm for CubitVector
408 [ + - ]: 68193 : double ai0 = N0 % T0;
409 [ + - ]: 68193 : double ai3 = N3 % T3;
410 : 68193 : double denom = 4 - ai * ai;
411 [ - + ]: 68193 : if( fabs( denom ) < 1e-20 )
412 : : {
413 : 0 : return MB_FAILURE; // CUBIT_FAILURE;
414 : : }
415 : 68193 : double row = 6.0e0 * ( 2.0e0 * ai0 + ai * ai3 ) / denom;
416 : 68193 : double omega = 6.0e0 * ( 2.0e0 * ai3 + ai * ai0 ) / denom;
417 [ + - ][ + - ]: 68193 : Vi[1] = Vi[0] + ( di * ( ( ( 6.0e0 * T0 ) - ( ( 2.0e0 * row ) * N0 ) + ( omega * N3 ) ) / 18.0e0 ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
418 [ + - ][ + - ]: 68193 : Vi[2] = Vi[3] - ( di * ( ( ( 6.0e0 * T3 ) + ( row * N0 ) - ( ( 2.0e0 * omega ) * N3 ) ) / 18.0e0 ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
419 : : // CartVect Wi[3];
420 : : // Wi[0] = Vi[1] - Vi[0];
421 : : // Wi[1] = Vi[2] - Vi[1];
422 : : // Wi[2] = Vi[3] - Vi[2];
423 : :
424 [ + - ][ + - ]: 68193 : Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
[ + - ]
425 [ + - ][ + - ]: 68193 : Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
[ + - ]
426 [ + - ][ + - ]: 68193 : Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
[ + - ]
427 : :
428 : 68193 : return MB_SUCCESS;
429 : : }
430 : :
431 : 44573 : ErrorCode SmoothFace::find_edges_orientations( EntityHandle edges[3], const EntityHandle* conn3,
432 : : int orient[3] ) // maybe we will set it?
433 : : {
434 : : // find the edge that is adjacent to 2 vertices at a time
435 [ + + ]: 178292 : for( int i = 0; i < 3; i++ )
436 : : {
437 : : // edge 0 is 1-2, 1 is 3-1, 2 is 0-1
438 : : EntityHandle v[2];
439 : 133719 : v[0] = conn3[( i + 1 ) % 3];
440 : 133719 : v[1] = conn3[( i + 2 ) % 3];
441 [ + - ]: 133719 : std::vector< EntityHandle > adjacencies;
442 : : // generate all edges for these two hexes
443 [ + - ]: 133719 : ErrorCode rval = _mb->get_adjacencies( v, 2, 1, false, adjacencies, Interface::INTERSECT );
444 [ - + ]: 133719 : if( MB_SUCCESS != rval ) return rval;
445 : :
446 : : // find the edge connected to both vertices, and then see its orientation
447 [ - + ]: 133719 : assert( adjacencies.size() == 1 );
448 : 133719 : const EntityHandle* conn2 = NULL;
449 : 133719 : int nnodes = 0;
450 [ + - ][ + - ]: 133719 : rval = _mb->get_connectivity( adjacencies[0], conn2, nnodes );
451 [ - + ]: 133719 : assert( rval == MB_SUCCESS );
452 [ - + ]: 133719 : assert( 2 == nnodes );
453 [ + - ]: 133719 : edges[i] = adjacencies[0];
454 : : // what is the story morning glory?
455 [ + + ][ + - ]: 133719 : if( conn2[0] == v[0] && conn2[1] == v[1] )
456 : 68177 : orient[i] = 1;
457 [ + - ][ + - ]: 65542 : else if( conn2[0] == v[1] && conn2[1] == v[0] )
458 : 65542 : orient[i] = -1;
459 : : else
460 [ + - ]: 133719 : return MB_FAILURE;
461 : 133719 : }
462 : 44573 : return MB_SUCCESS;
463 : : }
464 : 12 : ErrorCode SmoothFace::compute_internal_control_points_on_facets( double, Tag facetCtrlTag, Tag facetEdgeCtrlTag )
465 : : {
466 : : // collect from each triangle the control points in order
467 : : //
468 : :
469 : 12 : _facetCtrlTag = facetCtrlTag;
470 : 12 : _facetEdgeCtrlTag = facetEdgeCtrlTag;
471 : :
472 [ + - ][ + - ]: 44585 : for( Range::iterator it = _triangles.begin(); it != _triangles.end(); ++it )
[ + - ][ + - ]
[ + + ]
473 : : {
474 [ + - ]: 44573 : EntityHandle tri = *it;
475 : : // first get connectivity, and the edges
476 : : // we need a fast method to retrieve the adjacent edges to each triangle
477 : : const EntityHandle* conn3;
478 : : int nnodes;
479 [ + - ]: 44573 : ErrorCode rval = _mb->get_connectivity( tri, conn3, nnodes );
480 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
481 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
482 [ - + ]: 44573 : assert( 3 == nnodes );
483 : :
484 : : // would it be easier to do
485 [ + - ][ + + ]: 178292 : CartVect vNode[3]; // position at nodes
486 [ + - ]: 44573 : rval = _mb->get_coords( conn3, 3, (double*)&vNode[0] );
487 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
488 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
489 : :
490 : : // get gradients (normal) at each node of triangle
491 [ + - ][ + + ]: 178292 : CartVect NN[3];
492 [ + - ]: 44573 : rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
493 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
494 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
495 : :
496 : : EntityHandle edges[3];
497 : : int orient[3]; // + 1 or -1, if the edge is positive or negative within the face
498 [ + - ]: 44573 : rval = find_edges_orientations( edges, conn3, orient ); // maybe we will set it?
499 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
500 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
501 : : // maybe we will store some tags with edges and their orientation with respect to
502 : : // a triangle;
503 [ + - ][ + + ]: 846887 : CartVect P[3][5];
[ + + ]
504 [ + - ][ + + ]: 579449 : CartVect N[6], G[6];
[ + - ][ + + ]
505 : : // create the linear array for control points on edges, for storage (expensive !!!)
506 [ + - ][ + + ]: 445730 : CartVect CP[9];
507 : 44573 : int index = 0;
508 : : // maybe store a tag / entity handle for edges?
509 [ + + ]: 178292 : for( int i = 0; i < 3; i++ )
510 : : {
511 : : // populate P and N with the right vectors
512 : 133719 : int i1 = ( i + 1 ) % 3; // the first node of the edge
513 : 133719 : int i2 = ( i + 2 ) % 3; // the second node of the edge
514 : 133719 : N[2 * i] = NN[i1];
515 : 133719 : N[2 * i + 1] = NN[i2];
516 : 133719 : P[i][0] = vNode[i1];
517 [ + - ]: 133719 : rval = _mb->tag_get_data( _edgeCtrlTag, &edges[i], 1, &( P[i][1] ) );
518 : : // if sense is -1, swap 1 and 3 control points
519 [ + + ]: 133719 : if( orient[i] == -1 )
520 : : {
521 [ + - ]: 65542 : CartVect tmp;
522 : 65542 : tmp = P[i][1];
523 : 65542 : P[i][1] = P[i][3];
524 : 65542 : P[i][3] = tmp;
525 : : }
526 : 133719 : P[i][4] = vNode[i2];
527 [ + + ]: 534876 : for( int j = 1; j < 4; j++ )
528 : 401157 : CP[index++] = P[i][j];
529 : :
530 : : // the first edge control points
531 : : }
532 : :
533 : : // stat = facet->get_edge_control_points( P );
534 [ + - ]: 44573 : init_facet_control_points( N, P, G );
535 : : // what do we need to store in the tag control points?
536 [ + - ]: 44573 : rval = _mb->tag_set_data( _facetCtrlTag, &tri, 1, &G[0] );
537 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
538 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
539 : :
540 : : // store here again the 9 control points on the edges
541 [ + - ]: 44573 : rval = _mb->tag_set_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
542 [ - + ]: 44573 : assert( MB_SUCCESS == rval );
543 [ - + ]: 44573 : if( MB_SUCCESS != rval ) return rval;
544 : : // look at what we retrieve later
545 : :
546 : : // adjust the bounding box
547 : 44573 : int j = 0;
548 [ + + ]: 178292 : for( j = 0; j < 3; j++ )
549 [ + - ]: 133719 : adjust_bounding_box( vNode[j] );
550 : : // edge control points
551 [ + + ]: 445730 : for( j = 0; j < 9; j++ )
552 [ + - ]: 401157 : adjust_bounding_box( CP[j] );
553 : : // internal facet control points
554 [ + + ]: 312011 : for( j = 0; j < 6; j++ )
555 [ + - ]: 267438 : adjust_bounding_box( G[j] );
556 : : }
557 : 12 : return MB_SUCCESS;
558 : : }
559 : 802314 : void SmoothFace::adjust_bounding_box( CartVect& vect )
560 : : {
561 : : // _minim, _maxim
562 [ + + ]: 3209256 : for( int j = 0; j < 3; j++ )
563 : : {
564 [ + + ]: 2406942 : if( _minim[j] > vect[j] ) _minim[j] = vect[j];
565 [ + + ]: 2406942 : if( _maxim[j] < vect[j] ) _maxim[j] = vect[j];
566 : : }
567 : 802314 : }
568 : : //===============================================================
569 : : ////Function Name: init_facet_control_points
570 : : ////
571 : : ////Member Type: PRIVATE
572 : : ////Description: compute the control points for a facet
573 : : ////===============================================================
574 : 44573 : ErrorCode SmoothFace::init_facet_control_points( CartVect N[6], // vertex normals (per edge)
575 : : CartVect P[3][5], // edge control points
576 : : CartVect G[6] ) // return internal control points
577 : : {
578 [ + - ][ + + ]: 668595 : CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ][ + + ]
579 : : double denom;
580 : : double lambda[2], mu[2];
581 : :
582 : 44573 : ErrorCode rval = MB_SUCCESS;
583 : :
584 [ + + ]: 178292 : for( int i = 0; i < 3; i++ )
585 : : {
586 : 133719 : N0 = N[i * 2];
587 : 133719 : N3 = N[i * 2 + 1];
588 : 133719 : Vi[0] = P[i][0];
589 [ + - ][ + - ]: 133719 : Vi[1] = ( P[i][1] - 0.25 * P[i][0] ) / 0.75;
[ + - ]
590 [ + - ][ + - ]: 133719 : Vi[2] = ( P[i][3] - 0.25 * P[i][4] ) / 0.75;
[ + - ]
591 : 133719 : Vi[3] = P[i][4];
592 [ + - ]: 133719 : Wi[0] = Vi[1] - Vi[0];
593 [ + - ]: 133719 : Wi[1] = Vi[2] - Vi[1];
594 [ + - ]: 133719 : Wi[2] = Vi[3] - Vi[2];
595 [ + - ][ + - ]: 133719 : Di[0] = P[( i + 2 ) % 3][3] - 0.5 * ( P[i][1] + P[i][0] );
[ + - ]
596 [ + - ][ + - ]: 133719 : Di[3] = P[( i + 1 ) % 3][1] - 0.5 * ( P[i][4] + P[i][3] );
[ + - ]
597 [ + - ][ + - ]: 133719 : Ai[0] = ( N0 * Wi[0] ) / Wi[0].length();
[ + - ]
598 [ + - ][ + - ]: 133719 : Ai[2] = ( N3 * Wi[2] ) / Wi[2].length();
[ + - ]
599 [ + - ]: 133719 : Ai[1] = Ai[0] + Ai[2];
600 [ + - ]: 133719 : denom = Ai[1].length();
601 [ + - ]: 133719 : Ai[1] /= denom;
602 [ + - ][ + - ]: 133719 : lambda[0] = ( Di[0] % Wi[0] ) / ( Wi[0] % Wi[0] );
603 [ + - ][ + - ]: 133719 : lambda[1] = ( Di[3] % Wi[2] ) / ( Wi[2] % Wi[2] );
604 [ + - ]: 133719 : mu[0] = ( Di[0] % Ai[0] );
605 [ + - ]: 133719 : mu[1] = ( Di[3] % Ai[2] );
606 [ + - ][ + - ]: 133719 : G[i * 2] = 0.5 * ( P[i][1] + P[i][2] ) + 0.66666666666666 * lambda[0] * Wi[1] +
[ + - ][ + - ]
[ + - ]
607 [ + - ][ + - ]: 267438 : 0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0] * Ai[1] +
[ + - ]
608 [ + - ][ + - ]: 267438 : 0.33333333333333 * mu[1] * Ai[0];
609 [ + - ][ + - ]: 133719 : G[i * 2 + 1] = 0.5 * ( P[i][2] + P[i][3] ) + 0.33333333333333 * lambda[0] * Wi[2] +
[ + - ][ + - ]
[ + - ]
610 [ + - ][ + - ]: 267438 : 0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333 * mu[0] * Ai[2] +
[ + - ]
611 [ + - ][ + - ]: 267438 : 0.66666666666666 * mu[1] * Ai[1];
612 : : }
613 : 44573 : return rval;
614 : : }
615 : :
616 : 0 : void SmoothFace::DumpModelControlPoints()
617 : : {
618 : : // here, we will dump all control points from edges and facets (6 control points for each facet)
619 : : // we may also create some edges; maybe later...
620 : : // create a point3D file
621 : : // output a Point3D file (special visit format)
622 [ # # ]: 0 : unsigned long setId = _mb->id_from_handle( _set );
623 : 0 : char name[50] = { 0 };
624 : 0 : sprintf( name, "%lucontrol.Point3D", setId ); // name should be something 2control.Point3D
625 [ # # ]: 0 : std::ofstream point3DFile;
626 [ # # ][ # # ]: 0 : point3DFile.open( name ); //("control.Point3D");
627 [ # # ]: 0 : point3DFile << "# x y z \n";
628 [ # # ]: 0 : std::ofstream point3DEdgeFile;
629 : 0 : sprintf( name, "%lucontrolEdge.Point3D", setId ); //
630 [ # # ][ # # ]: 0 : point3DEdgeFile.open( name ); //("controlEdge.Point3D");
631 [ # # ]: 0 : point3DEdgeFile << "# x y z \n";
632 [ # # ]: 0 : std::ofstream smoothPoints;
633 : 0 : sprintf( name, "%lusmooth.Point3D", setId ); //
634 [ # # ][ # # ]: 0 : smoothPoints.open( name ); //("smooth.Point3D");
635 [ # # ]: 0 : smoothPoints << "# x y z \n";
636 [ # # ][ # # ]: 0 : CartVect controlPoints[3]; // edge control points
637 [ # # ][ # # ]: 0 : for( Range::iterator it = _edges.begin(); it != _edges.end(); ++it )
[ # # ][ # # ]
[ # # ]
638 : : {
639 [ # # ]: 0 : EntityHandle edge = *it;
640 : :
641 [ # # ]: 0 : _mb->tag_get_data( _edgeCtrlTag, &edge, 1, (double*)&controlPoints[0] );
642 [ # # ]: 0 : for( int i = 0; i < 3; i++ )
643 : : {
644 : 0 : CartVect& c = controlPoints[i];
645 [ # # ][ # # ]: 0 : point3DEdgeFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
646 : : }
647 : : }
648 [ # # ][ # # ]: 0 : CartVect controlTriPoints[6]; // triangle control points
649 [ # # ][ # # ]: 0 : CartVect P_facet[3]; // result in 3 "mid" control points
650 [ # # ][ # # ]: 0 : for( Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); ++it2 )
[ # # ][ # # ]
[ # # ]
651 : : {
652 [ # # ]: 0 : EntityHandle tri = *it2;
653 : :
654 [ # # ]: 0 : _mb->tag_get_data( _facetCtrlTag, &tri, 1, (double*)&controlTriPoints[0] );
655 : :
656 : : // draw a line of points between pairs of control points
657 : 0 : int numPoints = 7;
658 [ # # ]: 0 : for( int n = 0; n < numPoints; n++ )
659 : : {
660 : 0 : double a = 1. * n / ( numPoints - 1 );
661 : 0 : double b = 1.0 - a;
662 : :
663 [ # # ][ # # ]: 0 : P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
[ # # ]
664 : : // 1,2,1
665 [ # # ][ # # ]: 0 : P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
[ # # ]
666 : : // 1,1,2
667 [ # # ][ # # ]: 0 : P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
[ # # ]
668 [ # # ]: 0 : for( int i = 0; i < 3; i++ )
669 : : {
670 : 0 : CartVect& c = P_facet[i];
671 [ # # ][ # # ]: 0 : point3DFile << std::setprecision( 11 ) << c[0] << " " << c[1] << " " << c[2] << " \n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
672 : : }
673 : : }
674 : :
675 : : // evaluate for each triangle a lattice of points
676 : 0 : int N = 40;
677 [ # # ]: 0 : for( int k = 0; k <= N; k++ )
678 : : {
679 [ # # ]: 0 : for( int m = 0; m <= N - k; m++ )
680 : : {
681 : 0 : int n = N - m - k;
682 [ # # ]: 0 : CartVect areacoord( 1. * k / N, 1. * m / N, 1. * n / N );
683 [ # # ]: 0 : CartVect pt;
684 [ # # ]: 0 : eval_bezier_patch( tri, areacoord, pt );
685 [ # # ][ # # ]: 0 : smoothPoints << std::setprecision( 11 ) << pt[0] << " " << pt[1] << " " << pt[2] << " \n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
686 : : }
687 : : }
688 : : }
689 [ # # ]: 0 : point3DFile.close();
690 [ # # ]: 0 : smoothPoints.close();
691 [ # # ]: 0 : point3DEdgeFile.close();
692 : 0 : return;
693 : : }
694 : : //===========================================================================
695 : : // Function Name: evaluate_single
696 : : //
697 : : // Member Type: PUBLIC
698 : : // Description: evaluate edge not associated with a facet (this is used
699 : : // by camal edge mesher!!!)
700 : : // Note: t is a value from 0 to 1, for us
701 : : //===========================================================================
702 : 0 : ErrorCode SmoothFace::evaluate_smooth_edge( EntityHandle eh, double& tt, CartVect& outv )
703 : : {
704 [ # # ][ # # ]: 0 : CartVect P[2]; // P0 and P1
705 [ # # ][ # # ]: 0 : CartVect controlPoints[3]; // edge control points
706 : : double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
707 : :
708 : : // project the position to the linear edge
709 : : // t is from 0 to 1 only!!
710 : : // double tt = (t + 1) * 0.5;
711 [ # # ]: 0 : if( tt <= 0.0 ) tt = 0.0;
712 [ # # ]: 0 : if( tt >= 1.0 ) tt = 1.0;
713 : :
714 : 0 : int nnodes = 0;
715 : 0 : const EntityHandle* conn2 = NULL;
716 [ # # ]: 0 : ErrorCode rval = _mb->get_connectivity( eh, conn2, nnodes );
717 [ # # ]: 0 : assert( rval == MB_SUCCESS );
718 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
719 : :
720 [ # # ]: 0 : rval = _mb->get_coords( conn2, 2, (double*)&P[0] );
721 [ # # ]: 0 : assert( rval == MB_SUCCESS );
722 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
723 : :
724 [ # # ]: 0 : rval = _mb->tag_get_data( _edgeCtrlTag, &eh, 1, (double*)&controlPoints[0] );
725 [ # # ]: 0 : assert( rval == MB_SUCCESS );
726 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
727 : :
728 : 0 : t2 = tt * tt;
729 : 0 : t3 = t2 * tt;
730 : 0 : t4 = t3 * tt;
731 : 0 : one_minus_t = 1. - tt;
732 : 0 : one_minus_t2 = one_minus_t * one_minus_t;
733 : 0 : one_minus_t3 = one_minus_t2 * one_minus_t;
734 : 0 : one_minus_t4 = one_minus_t3 * one_minus_t;
735 : :
736 [ # # ][ # # ]: 0 : outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6. * one_minus_t2 * t2 * controlPoints[1] +
[ # # ][ # # ]
[ # # ][ # # ]
737 [ # # ][ # # ]: 0 : 4. * one_minus_t * t3 * controlPoints[2] + t4 * P[1];
[ # # ]
738 : :
739 : 0 : return MB_SUCCESS;
740 : : }
741 : :
742 : 37090 : ErrorCode SmoothFace::eval_bezier_patch( EntityHandle tri, CartVect& areacoord, CartVect& pt )
743 : : {
744 : : //
745 : : // interpolate internal control points
746 : :
747 [ + - ][ + + ]: 259630 : CartVect gctrl_pts[6];
748 : : // get the control points facet->get_control_points( gctrl_pts );
749 : : // init_facet_control_points( N, P, G) ;
750 : : // what do we need to store in the tag control points?
751 [ + - ]: 37090 : ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &tri, 1, &gctrl_pts[0] ); // get all 6 control points
752 [ - + ]: 37090 : assert( MB_SUCCESS == rval );
753 [ - + ]: 37090 : if( MB_SUCCESS != rval ) return rval;
754 : 37090 : const EntityHandle* conn3 = NULL;
755 : 37090 : int nnodes = 0;
756 [ + - ]: 37090 : rval = _mb->get_connectivity( tri, conn3, nnodes );
757 [ - + ]: 37090 : assert( MB_SUCCESS == rval );
758 : :
759 [ + - ][ + + ]: 148360 : CartVect vN[3];
760 [ + - ]: 37090 : _mb->get_coords( conn3, 3, (double*)&vN[0] ); // fill the coordinates of the vertices
761 : :
762 [ + - ][ + - ]: 37090 : if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
[ + + ]
763 : : {
764 : 2 : pt = vN[0];
765 : 2 : return MB_SUCCESS;
766 : : }
767 [ + - ][ + - ]: 37088 : if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
[ - + ]
768 : : {
769 : 0 : pt = vN[0];
770 : 0 : return MB_SUCCESS;
771 : : }
772 [ + - ][ + - ]: 37088 : if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
[ + + ]
773 : : {
774 : 4 : pt = vN[0];
775 : 4 : return MB_SUCCESS;
776 : : }
777 : :
778 [ + - ][ + + ]: 148336 : CartVect P_facet[3];
779 : : // 2,1,1
780 : : P_facet[0] =
781 [ + - ][ + - ]: 37084 : ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
782 : : // 1,2,1
783 : : P_facet[1] =
784 [ + - ][ + - ]: 37084 : ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
785 : : // 1,1,2
786 : : P_facet[2] =
787 [ + - ][ + - ]: 37084 : ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
788 : :
789 : : // sum the contribution from each of the control points
790 : :
791 [ + - ]: 37084 : pt = CartVect( 0. ); // set all to 0, we start adding / accumulating different parts
792 : : // first edge is from node 0 to 1, index 2 in
793 : :
794 : : // retrieve the points, in order, and the control points on edges
795 : :
796 : : // store here again the 9 control points on the edges
797 [ + - ][ + + ]: 370840 : CartVect CP[9];
798 [ + - ]: 37084 : rval = _mb->tag_get_data( _facetEdgeCtrlTag, &tri, 1, &CP[0] );
799 [ - + ]: 37084 : assert( MB_SUCCESS == rval );
800 : :
801 : : // CubitFacetEdge *edge;
802 : : // edge = facet->edge(2);! start with edge 2, from 0-1
803 : 37084 : int k = 0;
804 [ + - ][ + + ]: 222504 : CartVect ctrl_pts[5];
805 : : // edge->control_points(facet, ctrl_pts);
806 : 37084 : ctrl_pts[0] = vN[0]; //
807 [ + + ]: 148336 : for( k = 1; k < 4; k++ )
808 : 111252 : ctrl_pts[k] = CP[k + 5]; // for edge index 2
809 : 37084 : ctrl_pts[4] = vN[1]; //
810 : :
811 : : // i=4; j=0; k=0;
812 [ + - ][ + - ]: 37084 : double B = mbquart( areacoord[0] );
[ + - ][ + - ]
813 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[0];
814 : :
815 : : // i=3; j=1; k=0;
816 [ + - ][ + - ]: 37084 : B = 4.0 * mbcube( areacoord[0] ) * areacoord[1];
[ + - ][ + - ]
817 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[1];
818 : :
819 : : // i=2; j=2; k=0;
820 [ + - ][ + - ]: 37084 : B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[1] );
[ + - ][ + - ]
821 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[2];
822 : :
823 : : // i=1; j=3; k=0;
824 [ + - ][ + - ]: 37084 : B = 4.0 * areacoord[0] * mbcube( areacoord[1] );
[ + - ][ + - ]
825 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[3];
826 : :
827 : : // edge = facet->edge(0);
828 : : // edge->control_points(facet, ctrl_pts);
829 : : // edge index 0, from 1 to 2
830 : 37084 : ctrl_pts[0] = vN[1]; //
831 [ + + ]: 148336 : for( k = 1; k < 4; k++ )
832 : 111252 : ctrl_pts[k] = CP[k - 1]; // for edge index 0
833 : 37084 : ctrl_pts[4] = vN[2]; //
834 : :
835 : : // i=0; j=4; k=0;
836 [ + - ][ + - ]: 37084 : B = mbquart( areacoord[1] );
[ + - ][ + - ]
837 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[0];
838 : :
839 : : // i=0; j=3; k=1;
840 [ + - ][ + - ]: 37084 : B = 4.0 * mbcube( areacoord[1] ) * areacoord[2];
[ + - ][ + - ]
841 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[1];
842 : :
843 : : // i=0; j=2; k=2;
844 [ + - ][ + - ]: 37084 : B = 6.0 * mbsqr( areacoord[1] ) * mbsqr( areacoord[2] );
[ + - ][ + - ]
845 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[2];
846 : :
847 : : // i=0; j=1; k=3;
848 [ + - ][ + - ]: 37084 : B = 4.0 * areacoord[1] * mbcube( areacoord[2] );
[ + - ][ + - ]
849 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[3];
850 : :
851 : : // edge = facet->edge(1);
852 : : // edge->control_points(facet, ctrl_pts);
853 : : // edge index 1, from 2 to 0
854 : 37084 : ctrl_pts[0] = vN[2]; //
855 [ + + ]: 148336 : for( k = 1; k < 4; k++ )
856 : 111252 : ctrl_pts[k] = CP[k + 2]; // for edge index 0
857 : 37084 : ctrl_pts[4] = vN[0]; //
858 : :
859 : : // i=0; j=0; k=4;
860 [ + - ][ + - ]: 37084 : B = mbquart( areacoord[2] );
[ + - ][ + - ]
861 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[0];
862 : :
863 : : // i=1; j=0; k=3;
864 [ + - ][ + - ]: 37084 : B = 4.0 * areacoord[0] * mbcube( areacoord[2] );
[ + - ][ + - ]
865 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[1];
866 : :
867 : : // i=2; j=0; k=2;
868 [ + - ][ + - ]: 37084 : B = 6.0 * mbsqr( areacoord[0] ) * mbsqr( areacoord[2] );
[ + - ][ + - ]
869 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[2];
870 : :
871 : : // i=3; j=0; k=1;
872 [ + - ][ + - ]: 37084 : B = 4.0 * mbcube( areacoord[0] ) * areacoord[2];
[ + - ][ + - ]
873 [ + - ][ + - ]: 37084 : pt += B * ctrl_pts[3];
874 : :
875 : : // i=2; j=1; k=1;
876 [ + - ][ + - ]: 37084 : B = 12.0 * mbsqr( areacoord[0] ) * areacoord[1] * areacoord[2];
[ + - ][ + - ]
877 [ + - ][ + - ]: 37084 : pt += B * P_facet[0];
878 : :
879 : : // i=1; j=2; k=1;
880 [ + - ][ + - ]: 37084 : B = 12.0 * areacoord[0] * mbsqr( areacoord[1] ) * areacoord[2];
[ + - ][ + - ]
881 [ + - ][ + - ]: 37084 : pt += B * P_facet[1];
882 : :
883 : : // i=1; j=1; k=2;
884 [ + - ][ + - ]: 37084 : B = 12.0 * areacoord[0] * areacoord[1] * mbsqr( areacoord[2] );
[ + - ][ + - ]
885 [ + - ][ + - ]: 37084 : pt += B * P_facet[2];
886 : :
887 : 37090 : return MB_SUCCESS;
888 : : }
889 : :
890 : : //===========================================================================
891 : : // Function Name: project_to_facet_plane
892 : : //
893 : : // Member Type: PUBLIC
894 : : // Descriptoin: Project a point to the plane of a facet
895 : : //===========================================================================
896 : 2680 : void SmoothFace::project_to_facet_plane( EntityHandle tri, CartVect& pt, CartVect& point_on_plane,
897 : : double& dist_to_plane )
898 : : {
899 : : double plane[4];
900 : :
901 [ + - ]: 2680 : ErrorCode rval = _mb->tag_get_data( _planeTag, &tri, 1, plane );
902 [ - + ]: 2680 : if( MB_SUCCESS != rval ) return;
903 [ - + ]: 2680 : assert( rval == MB_SUCCESS );
904 : : // _planeTag
905 [ + - ]: 2680 : CartVect normal( &plane[0] ); // just first 3 components are used
906 : :
907 [ + - ]: 2680 : double dist = normal % pt + plane[3]; // coeff d is saved!!!
908 : 2680 : dist_to_plane = fabs( dist );
909 [ + - ][ + - ]: 2680 : point_on_plane = pt - dist * normal;
910 : :
911 : 2680 : return;
912 : : }
913 : :
914 : : //===========================================================================
915 : : // Function Name: facet_area_coordinate
916 : : //
917 : : // Member Type: PUBLIC
918 : : // Descriptoin: Determine the area coordinates of a point on the plane
919 : : // of a facet
920 : : //===========================================================================
921 : 2680 : void SmoothFace::facet_area_coordinate( EntityHandle facet, CartVect& pt_on_plane, CartVect& areacoord )
922 : : {
923 : :
924 : 2680 : const EntityHandle* conn3 = NULL;
925 : 2680 : int nnodes = 0;
926 [ + - ]: 2680 : ErrorCode rval = _mb->get_connectivity( facet, conn3, nnodes );
927 [ - + ]: 2680 : assert( MB_SUCCESS == rval );
928 : : if( rval ) {} // empty statement to prevent compiler warning
929 : :
930 : : // double coords[9]; // store the coordinates for the nodes
931 : : //_mb->get_coords(conn3, 3, coords);
932 [ + - ][ + + ]: 10720 : CartVect p[3];
933 [ + - ]: 2680 : rval = _mb->get_coords( conn3, 3, (double*)&p[0] );
934 [ - + ]: 2680 : assert( MB_SUCCESS == rval );
935 : : if( rval ) {} // empty statement to prevent compiler warning
936 : : double plane[4];
937 : :
938 [ + - ]: 2680 : rval = _mb->tag_get_data( _planeTag, &facet, 1, plane );
939 [ - + ]: 2680 : assert( rval == MB_SUCCESS );
940 : : if( rval ) {} // empty statement to prevent compiler warning
941 [ + - ]: 2680 : CartVect normal( &plane[0] ); // just first 3 components are used
942 : :
943 : : double area2;
944 : :
945 : 2680 : double tol = GEOMETRY_RESABS * 1.e-5; // 1.e-11;
946 : :
947 [ + - ]: 2680 : CartVect v1( p[1] - p[0] );
948 [ + - ]: 2680 : CartVect v2( p[2] - p[0] );
949 : :
950 [ + - ][ + - ]: 2680 : area2 = ( v1 * v2 ).length_squared(); // the same for CartVect
951 [ - + ]: 2680 : if( area2 < 100 * tol ) { tol = .01 * area2; }
952 [ + - ][ + - ]: 2680 : CartVect absnorm( fabs( normal[0] ), fabs( normal[1] ), fabs( normal[2] ) );
[ + - ][ + - ]
953 : :
954 : : // project to the closest coordinate plane so we only have to do this in 2D
955 : :
956 [ + - ][ + - ]: 2680 : if( absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2] )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
957 : : {
958 [ + - ][ + - ]: 85 : area2 = determ3( p[0][1], p[0][2], p[1][1], p[1][2], p[2][1], p[2][2] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
959 [ - + ]: 85 : if( fabs( area2 ) < tol )
960 : : {
961 [ # # ]: 0 : areacoord = CartVect( -std::numeric_limits< double >::min() ); // .set(
962 : : // -std::numeric_limits<double>::min(),
963 : : // -std::numeric_limits<double>::min(),
964 : : // -std::numeric_limits<double>::min() );
965 : : }
966 [ + - ][ - + ]: 85 : else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
967 : : {
968 [ # # ]: 0 : areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
969 : : }
970 [ + - ][ - + ]: 85 : else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
971 : : {
972 [ # # ]: 0 : areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
973 : : }
974 [ + - ][ - + ]: 85 : else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
975 : : {
976 [ # # ]: 0 : areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
977 : : }
978 : : else
979 : : {
980 : :
981 [ + - ][ + - ]: 85 : areacoord[0] = determ3( pt_on_plane[1], pt_on_plane[2], p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
982 : :
983 [ + - ][ + - ]: 85 : areacoord[1] = determ3( p[0][1], p[0][2], pt_on_plane[1], pt_on_plane[2], p[2][1], p[2][2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
984 : :
985 [ + - ][ + - ]: 85 : areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2], pt_on_plane[1], pt_on_plane[2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
986 : : }
987 : : }
988 [ + - ][ + - ]: 2595 : else if( absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2] )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
989 : : {
990 : :
991 [ + - ][ + - ]: 89 : area2 = determ3( p[0][0], p[0][2], p[1][0], p[1][2], p[2][0], p[2][2] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
992 [ - + ]: 89 : if( fabs( area2 ) < tol )
993 : : {
994 [ # # ]: 0 : areacoord = CartVect( -std::numeric_limits< double >::min() ); //.set(
995 : : //-std::numeric_limits<double>::min(),
996 : : //-std::numeric_limits<double>::min(),
997 : : //-std::numeric_limits<double>::min() );
998 : : }
999 [ + - ][ - + ]: 89 : else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
1000 : : {
1001 [ # # ]: 0 : areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
1002 : : }
1003 [ + - ][ - + ]: 89 : else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
1004 : : {
1005 [ # # ]: 0 : areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
1006 : : }
1007 [ + - ][ - + ]: 89 : else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
1008 : : {
1009 [ # # ]: 0 : areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
1010 : : }
1011 : : else
1012 : : {
1013 : :
1014 [ + - ][ + - ]: 89 : areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[2], p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1015 : :
1016 [ + - ][ + - ]: 89 : areacoord[1] = determ3( p[0][0], p[0][2], pt_on_plane[0], pt_on_plane[2], p[2][0], p[2][2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1017 : :
1018 [ + - ][ + - ]: 89 : areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2], pt_on_plane[0], pt_on_plane[2] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1019 : : }
1020 : : }
1021 : : else
1022 : : {
1023 : : /*area2 = determ3(pt0->x(), pt0->y(),
1024 : : pt1->x(), pt1->y(),
1025 : : pt2->x(), pt2->y());*/
1026 [ + - ][ + - ]: 2506 : area2 = determ3( p[0][0], p[0][1], p[1][0], p[1][1], p[2][0], p[2][1] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1027 [ - + ]: 2506 : if( fabs( area2 ) < tol )
1028 : : {
1029 [ # # ]: 0 : areacoord = CartVect( -std::numeric_limits< double >::min() ); //.set(
1030 : : //-std::numeric_limits<double>::min(),
1031 : : //-std::numeric_limits<double>::min(),
1032 : : //-std::numeric_limits<double>::min() );
1033 : : }
1034 [ + - ][ - + ]: 2506 : else if( within_tolerance( p[0], pt_on_plane, GEOMETRY_RESABS ) )
1035 : : {
1036 [ # # ]: 0 : areacoord = CartVect( 1., 0., 0. ); //.set( 1.0, 0.0, 0.0 );
1037 : : }
1038 [ + - ][ - + ]: 2506 : else if( within_tolerance( p[1], pt_on_plane, GEOMETRY_RESABS ) )
1039 : : {
1040 [ # # ]: 0 : areacoord = CartVect( 0., 1., 0. ); //.set( 0.0, 1.0, 0.0 );
1041 : : }
1042 [ + - ][ - + ]: 2506 : else if( within_tolerance( p[2], pt_on_plane, GEOMETRY_RESABS ) )
1043 : : {
1044 [ # # ]: 0 : areacoord = CartVect( 0., 0., 1. ); //.set( 0.0, 0.0, 1.0 );
1045 : : }
1046 : : else
1047 : : {
1048 : :
1049 [ + - ][ + - ]: 2506 : areacoord[0] = determ3( pt_on_plane[0], pt_on_plane[1], p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1050 : :
1051 [ + - ][ + - ]: 2506 : areacoord[1] = determ3( p[0][0], p[0][1], pt_on_plane[0], pt_on_plane[1], p[2][0], p[2][1] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1052 : :
1053 [ + - ][ + - ]: 2506 : areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1], pt_on_plane[0], pt_on_plane[1] ) / area2;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1054 : : }
1055 : : }
1056 : 2680 : }
1057 : :
1058 : 1777 : ErrorCode SmoothFace::project_to_facets_main( CartVect& this_point, bool trim, bool& outside,
1059 : : CartVect* closest_point_ptr, CartVect* normal_ptr )
1060 : : {
1061 : :
1062 : : // if there are a lot of facets on this surface - use the OBB search first
1063 : : // to narrow the selection
1064 : :
1065 : 1777 : _evaluationsCounter++;
1066 : 1777 : double tolerance = 1.e-5;
1067 [ + - ]: 1777 : std::vector< EntityHandle > facets_out;
1068 : : // we will start with a list of facets anyway, the best among them wins
1069 : : ErrorCode rval =
1070 [ + - ][ + - ]: 1777 : _my_geomTopoTool->obb_tree()->closest_to_location( (double*)&this_point, _obb_root, tolerance, facets_out );
1071 [ - + ]: 1777 : if( MB_SUCCESS != rval ) return rval;
1072 : :
1073 : 1777 : int interpOrder = 4;
1074 : 1777 : double compareTol = 1.e-5;
1075 [ + - ]: 1777 : EntityHandle lastFacet = facets_out.front();
1076 : : rval = project_to_facets( facets_out, lastFacet, interpOrder, compareTol, this_point, trim, outside,
1077 [ + - ]: 1777 : closest_point_ptr, normal_ptr );
1078 : :
1079 : 1777 : return rval;
1080 : : }
1081 : 1777 : ErrorCode SmoothFace::project_to_facets( std::vector< EntityHandle >& facet_list, EntityHandle& lastFacet,
1082 : : int interpOrder, double compareTol, CartVect& this_point, bool, bool& outside,
1083 : : CartVect* closest_point_ptr, CartVect* normal_ptr )
1084 : : {
1085 : :
1086 : 1777 : bool outside_facet = false;
1087 : 1777 : bool best_outside_facet = true;
1088 : 1777 : double mindist = 1.e20;
1089 [ + - ][ + - ]: 1777 : CartVect close_point, best_point( mindist, mindist, mindist ), best_areacoord;
[ + - ]
1090 : 1777 : EntityHandle best_facet = 0L; // no best facet found yet
1091 : : EntityHandle facet;
1092 [ - + ]: 1777 : assert( facet_list.size() > 0 );
1093 : :
1094 : 1777 : double big_dist = compareTol * 1.0e3;
1095 : :
1096 : : // from the list of close facets, determine the closest point
1097 [ + + ]: 4457 : for( size_t i = 0; i < facet_list.size(); i++ )
1098 : : {
1099 [ + - ]: 2680 : facet = facet_list[i];
1100 [ + - ]: 2680 : CartVect pt_on_plane;
1101 : : double dist_to_plane;
1102 [ + - ]: 2680 : project_to_facet_plane( facet, this_point, pt_on_plane, dist_to_plane );
1103 : :
1104 [ + - ]: 2680 : CartVect areacoord;
1105 : : // CartVect close_point;
1106 [ + - ]: 2680 : facet_area_coordinate( facet, pt_on_plane, areacoord );
1107 [ + - ]: 2680 : if( interpOrder != 0 )
1108 : : {
1109 : :
1110 : : // modify the areacoord - project to the bezier patch- snaps to the
1111 : : // edge of the patch if necessary
1112 : :
1113 [ + - ][ - + ]: 2680 : if( project_to_facet( facet, this_point, areacoord, close_point, outside_facet, compareTol ) != MB_SUCCESS )
1114 : 0 : { return MB_FAILURE; }
1115 : : // if (closest_point_ptr)
1116 : : //*closest_point_ptr = close_point;
1117 : : }
1118 : : // keep track of the minimum distance
1119 : :
1120 [ + - ][ + - ]: 2680 : double dist = ( close_point - this_point ).length(); // close_point.distance_between(this_point);
1121 [ + + ][ + + ]: 2680 : if( ( best_outside_facet == outside_facet && dist < mindist ) ||
[ + + ]
1122 [ + + ][ + + ]: 1568 : ( best_outside_facet && !outside_facet && ( dist < big_dist || best_facet == 0L /*!best_facet*/ ) ) )
[ + - ]
1123 : : {
1124 : 2236 : mindist = dist;
1125 : 2236 : best_point = close_point;
1126 : 2236 : best_facet = facet;
1127 : 2236 : best_areacoord = areacoord;
1128 : 2236 : best_outside_facet = outside_facet;
1129 : :
1130 [ - + ]: 2236 : if( dist < compareTol ) { break; }
1131 : 2236 : big_dist = 10.0 * mindist;
1132 : : }
1133 : : // facet->marked(1);
1134 : : // used_facet_list.append(facet);
1135 : : }
1136 : :
1137 [ + + ]: 1777 : if( normal_ptr )
1138 : : {
1139 [ + - ]: 14 : CartVect normal;
1140 [ + - ][ - + ]: 14 : if( eval_bezier_patch_normal( best_facet, best_areacoord, normal ) != MB_SUCCESS ) { return MB_FAILURE; }
1141 : 14 : *normal_ptr = normal;
1142 : : }
1143 : :
1144 [ + + ]: 1777 : if( closest_point_ptr ) { *closest_point_ptr = best_point; }
1145 : :
1146 : 1777 : outside = best_outside_facet;
1147 : 1777 : lastFacet = best_facet;
1148 : :
1149 : 1777 : return MB_SUCCESS;
1150 : : // end copy
1151 : : }
1152 : :
1153 : : //===========================================================================
1154 : : // Function Name: project_to_patch
1155 : : //
1156 : : // Member Type: PUBLIC
1157 : : // Description: Project a point to a bezier patch. Pass in the areacoord
1158 : : // of the point projected to the linear facet. Function
1159 : : // assumes that the point is contained within the patch -
1160 : : // if not, it will project to one of its edges.
1161 : : //===========================================================================
1162 : 2680 : ErrorCode SmoothFace::project_to_patch( EntityHandle facet, // (IN) the facet where the patch is defined
1163 : : CartVect& ac, // (IN) area coordinate initial guess (from linear facet)
1164 : : CartVect& pt, // (IN) point we are projecting to patch
1165 : : CartVect& eval_pt, // (OUT) The projected point
1166 : : CartVect* eval_norm, // (OUT) normal at evaluated point
1167 : : bool& outside, // (OUT) the closest point on patch to pt is on an edge
1168 : : double compare_tol, // (IN) comparison tolerance
1169 : : int edge_id ) // (IN) only used if this is to be projected to one
1170 : : // of the edges. Otherwise, should be -1
1171 : : {
1172 : 2680 : ErrorCode status = MB_SUCCESS;
1173 : :
1174 : : // see if we are at a vertex
1175 : :
1176 : : #define INCR 0.01
1177 : 2680 : const double tol = compare_tol;
1178 : :
1179 [ + - ][ - + ]: 2680 : if( is_at_vertex( facet, pt, ac, compare_tol, eval_pt, eval_norm ) )
1180 : : {
1181 : 0 : outside = false;
1182 : 0 : return MB_SUCCESS;
1183 : : }
1184 : :
1185 : : // check if the start ac is inside the patch -if not, then move it there
1186 : :
1187 : 2680 : int nout = 0;
1188 : 2680 : const double atol = 0.001;
1189 [ + - ][ + + ]: 2680 : if( move_ac_inside( ac, atol ) ) nout++;
1190 : :
1191 : 2680 : int diverge = 0;
1192 : 2680 : int iter = 0;
1193 [ + - ]: 2680 : CartVect newpt;
1194 [ + - ]: 2680 : eval_bezier_patch( facet, ac, newpt );
1195 [ + - ]: 2680 : CartVect move = pt - newpt;
1196 [ + - ]: 2680 : double lastdist = move.length();
1197 : 2680 : double bestdist = lastdist;
1198 : 2680 : CartVect bestac = ac;
1199 : 2680 : CartVect bestpt = newpt;
1200 [ + - ]: 2680 : CartVect bestnorm( 0, 0, 0 );
1201 : :
1202 : : // If we are already close enough, then return now
1203 : :
1204 [ - + ][ # # ]: 2680 : if( lastdist <= tol && !eval_norm && nout == 0 )
[ # # ]
1205 : : {
1206 : 0 : eval_pt = pt;
1207 : 0 : outside = false;
1208 : 0 : return status;
1209 : : }
1210 : :
1211 : : double ratio, mag, umove, vmove, det, distnew, movedist;
1212 : 2680 : CartVect lastpt = newpt;
1213 : 2680 : CartVect lastac = ac;
1214 [ + - ]: 2680 : CartVect norm;
1215 [ + - ][ + - ]: 2680 : CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1216 [ + - ][ + - ]: 2680 : CartVect du, dv, newac;
[ + - ]
1217 : 2680 : bool done = false;
1218 [ + + ]: 14150 : while( !done )
1219 : : {
1220 : :
1221 : : // We will be locating the projected point within the u,v,w coordinate
1222 : : // system of the triangular bezier patch. Since u+v+w=1, the components
1223 : : // are not linearly independent. We will choose only two of the
1224 : : // coordinates to use and compute the third.
1225 : :
1226 : : int system;
1227 [ + - ][ + - ]: 11470 : if( lastac[0] >= lastac[1] && lastac[0] >= lastac[2] ) { system = 0; }
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
1228 [ + - ][ + - ]: 7641 : else if( lastac[1] >= lastac[2] )
[ + + ]
1229 : : {
1230 : 4011 : system = 1;
1231 : : }
1232 : : else
1233 : : {
1234 : 3630 : system = 2;
1235 : : }
1236 : :
1237 : : // compute the surface derivatives with respect to each
1238 : : // of the barycentric coordinates
1239 : :
1240 [ + + ][ + + ]: 11470 : if( system == 1 || system == 2 )
1241 : : {
1242 [ + - ][ + - ]: 7641 : xac[0] = lastac[0] + INCR; // xac.x( lastac.x() + INCR );
1243 [ + - ][ + - ]: 7641 : if( lastac[1] + lastac[2] == 0.0 ) return MB_FAILURE;
[ - + ]
1244 [ + - ][ + - ]: 7641 : ratio = lastac[2] / ( lastac[1] + lastac[2] ); // ratio = lastac.z() / (lastac.y() + lastac.z());
[ + - ]
1245 [ + - ][ + - ]: 7641 : xac[1] = ( 1.0 - xac[0] ) * ( 1.0 - ratio ); // xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
1246 [ + - ][ + - ]: 7641 : xac[2] = 1.0 - xac[0] - xac[1]; // xac.z( 1.0 - xac.x() - xac.y() );
[ + - ]
1247 [ + - ]: 7641 : eval_bezier_patch( facet, xac, xpt );
1248 [ + - ]: 7641 : xvec = xpt - lastpt;
1249 [ + - ]: 7641 : xvec /= INCR;
1250 : : }
1251 [ + + ][ + + ]: 11470 : if( system == 0 || system == 2 )
1252 : : {
1253 [ + - ][ + - ]: 7459 : yac[1] = ( lastac[1] + INCR ); // yac.y( lastac.y() + INCR );
1254 [ + - ][ + - ]: 7459 : if( lastac[0] + lastac[2] == 0.0 ) // if (lastac.x() + lastac.z() == 0.0)
[ - + ]
1255 : 0 : return MB_FAILURE;
1256 [ + - ][ + - ]: 7459 : ratio = lastac[2] / ( lastac[0] + lastac[2] ); // ratio = lastac.z() / (lastac.x() + lastac.z());
[ + - ]
1257 [ + - ][ + - ]: 7459 : yac[0] = ( ( 1.0 - yac[1] ) * ( 1.0 - ratio ) ); // yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
1258 [ + - ][ + - ]: 7459 : yac[2] = ( 1.0 - yac[0] - yac[1] ); // yac.z( 1.0 - yac.x() - yac.y() );
[ + - ]
1259 [ + - ]: 7459 : eval_bezier_patch( facet, yac, ypt );
1260 [ + - ]: 7459 : yvec = ypt - lastpt;
1261 [ + - ]: 7459 : yvec /= INCR;
1262 : : }
1263 [ + + ][ + + ]: 11470 : if( system == 0 || system == 1 )
1264 : : {
1265 [ + - ][ + - ]: 7840 : zac[2] = ( lastac[2] + INCR ); // zac.z( lastac.z() + INCR );
1266 [ + - ][ + - ]: 7840 : if( lastac[0] + lastac[1] == 0.0 ) // if (lastac.x() + lastac.y() == 0.0)
[ - + ]
1267 : 0 : return MB_FAILURE;
1268 [ + - ][ + - ]: 7840 : ratio = lastac[1] / ( lastac[0] + lastac[1] ); // ratio = lastac.y() / (lastac.x() + lastac.y());
[ + - ]
1269 [ + - ][ + - ]: 7840 : zac[0] = ( ( 1.0 - zac[2] ) * ( 1.0 - ratio ) ); // zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
1270 [ + - ][ + - ]: 7840 : zac[1] = ( 1.0 - zac[0] - zac[2] ); // zac.y( 1.0 - zac.x() - zac.z() );
[ + - ]
1271 [ + - ]: 7840 : eval_bezier_patch( facet, zac, zpt );
1272 [ + - ]: 7840 : zvec = zpt - lastpt;
1273 [ + - ]: 7840 : zvec /= INCR;
1274 : : }
1275 : :
1276 : : // compute the surface normal
1277 : :
1278 [ + + + - ]: 11470 : switch( system )
1279 : : {
1280 : : case 0:
1281 : 3829 : du = yvec;
1282 : 3829 : dv = zvec;
1283 : 3829 : break;
1284 : : case 1:
1285 : 4011 : du = zvec;
1286 : 4011 : dv = xvec;
1287 : 4011 : break;
1288 : : case 2:
1289 : 3630 : du = xvec;
1290 : 3630 : dv = yvec;
1291 : 3630 : break;
1292 : : }
1293 [ + - ]: 11470 : norm = du * dv;
1294 [ + - ]: 11470 : mag = norm.length();
1295 [ - + ]: 11470 : if( mag < DBL_EPSILON )
1296 : : {
1297 : 0 : return MB_FAILURE;
1298 : : // do something else here (it is likely a flat triangle -
1299 : : // so try evaluating just an edge of the bezier patch)
1300 : : }
1301 [ + - ]: 11470 : norm /= mag;
1302 [ + + ]: 11470 : if( iter == 0 ) bestnorm = norm;
1303 : :
1304 : : // project the move vector to the tangent plane
1305 : :
1306 [ + - ][ + - ]: 11470 : move = ( norm * move ) * norm;
1307 : :
1308 : : // compute an equivalent u-v-w vector
1309 : :
1310 [ + - ][ + - ]: 11470 : CartVect absnorm( fabs( norm[0] ), fabs( norm[1] ), fabs( norm[2] ) );
[ + - ][ + - ]
1311 [ + - ][ + - ]: 11470 : if( absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0] )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
1312 : : {
1313 [ + - ][ + - ]: 10674 : det = du[0] * dv[1] - dv[0] * du[1];
[ + - ][ + - ]
1314 [ - + ]: 10674 : if( fabs( det ) <= DBL_EPSILON )
1315 : : {
1316 : 0 : return MB_FAILURE; // do something else here
1317 : : }
1318 [ + - ][ + - ]: 10674 : umove = ( move[0] * dv[1] - dv[0] * move[1] ) / det;
[ + - ][ + - ]
1319 [ + - ][ + - ]: 10674 : vmove = ( du[0] * move[1] - move[0] * du[1] ) / det;
[ + - ][ + - ]
1320 : : }
1321 [ + - ][ + - ]: 796 : else if( absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0] )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
1322 : : {
1323 [ + - ][ + - ]: 351 : det = du[0] * dv[2] - dv[0] * du[2];
[ + - ][ + - ]
1324 [ - + ]: 351 : if( fabs( det ) <= DBL_EPSILON ) { return MB_FAILURE; }
1325 [ + - ][ + - ]: 351 : umove = ( move[0] * dv[2] - dv[0] * move[2] ) / det;
[ + - ][ + - ]
1326 [ + - ][ + - ]: 351 : vmove = ( du[0] * move[2] - move[0] * du[2] ) / det;
[ + - ][ + - ]
1327 : : }
1328 : : else
1329 : : {
1330 [ + - ][ + - ]: 445 : det = du[1] * dv[2] - dv[1] * du[2];
[ + - ][ + - ]
1331 [ - + ]: 445 : if( fabs( det ) <= DBL_EPSILON ) { return MB_FAILURE; }
1332 [ + - ][ + - ]: 445 : umove = ( move[1] * dv[2] - dv[1] * move[2] ) / det;
[ + - ][ + - ]
1333 [ + - ][ + - ]: 445 : vmove = ( du[1] * move[2] - move[1] * du[2] ) / det;
[ + - ][ + - ]
1334 : : }
1335 : :
1336 : : /* === compute the new u-v coords and evaluate surface at new location */
1337 : :
1338 [ + + + - ]: 11470 : switch( system )
1339 : : {
1340 : : case 0:
1341 [ + - ][ + - ]: 3829 : newac[1] = ( lastac[1] + umove ); // newac.y( lastac.y() + umove );
1342 [ + - ][ + - ]: 3829 : newac[2] = ( lastac[2] + vmove ); // newac.z( lastac.z() + vmove );
1343 [ + - ][ + - ]: 3829 : newac[0] = ( 1.0 - newac[1] - newac[2] ); // newac.x( 1.0 - newac.y() - newac.z() );
[ + - ]
1344 : 3829 : break;
1345 : : case 1:
1346 [ + - ][ + - ]: 4011 : newac[2] = ( lastac[2] + umove ); // newac.z( lastac.z() + umove );
1347 [ + - ][ + - ]: 4011 : newac[0] = ( lastac[0] + vmove ); // newac.x( lastac.x() + vmove );
1348 [ + - ][ + - ]: 4011 : newac[1] = ( 1.0 - newac[2] - newac[0] ); // newac.y( 1.0 - newac.z() - newac.x() );
[ + - ]
1349 : 4011 : break;
1350 : : case 2:
1351 [ + - ][ + - ]: 3630 : newac[0] = ( lastac[0] + umove ); // newac.x( lastac.x() + umove );
1352 [ + - ][ + - ]: 3630 : newac[1] = ( lastac[1] + vmove ); // newac.y( lastac.y() + vmove );
1353 [ + - ][ + - ]: 3630 : newac[2] = ( 1.0 - newac[0] - newac[1] ); // newac.z( 1.0 - newac.x() - newac.y() );
[ + - ]
1354 : 3630 : break;
1355 : : }
1356 : :
1357 : : // Keep it inside the patch
1358 : :
1359 [ + - ][ + + ]: 11470 : if( newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol ) { nout = 0; }
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ]
1360 : : else
1361 : : {
1362 [ + - ][ + - ]: 2110 : if( move_ac_inside( newac, atol ) ) nout++;
1363 : : }
1364 : :
1365 : : // Evaluate at the new location
1366 : :
1367 [ - + ][ # # ]: 11470 : if( edge_id != -1 ) ac_at_edge( newac, newac, edge_id ); // move to edge first
1368 [ + - ]: 11470 : eval_bezier_patch( facet, newac, newpt );
1369 : :
1370 : : // Check for convergence
1371 : :
1372 [ + - ][ + - ]: 11470 : distnew = ( pt - newpt ).length(); // pt.distance_between(newpt);
1373 [ + - ]: 11470 : move = newpt - lastpt;
1374 [ + - ]: 11470 : movedist = move.length();
1375 [ + + ][ - + ]: 11470 : if( movedist < tol || distnew < tol )
1376 : : {
1377 : 2543 : done = true;
1378 [ + + ]: 3269 : if( distnew < bestdist )
1379 : : {
1380 : 726 : bestdist = distnew;
1381 : 726 : bestac = newac;
1382 : 726 : bestpt = newpt;
1383 : 726 : bestnorm = norm;
1384 : : }
1385 : : }
1386 : : else
1387 : : {
1388 : :
1389 : : // don't allow more than 30 iterations
1390 : :
1391 : 8927 : iter++;
1392 [ - + ]: 8927 : if( iter > 30 )
1393 : : {
1394 : : // if (movedist > tol * 100.0) nout=1;
1395 : 0 : done = true;
1396 : : }
1397 : :
1398 : : // Check for divergence - don't allow more than 5 divergent
1399 : : // iterations
1400 : :
1401 [ + + ]: 8927 : if( distnew > lastdist )
1402 : : {
1403 : 2177 : diverge++;
1404 [ + + ]: 2177 : if( diverge > 10 )
1405 : : {
1406 : 4 : done = true;
1407 : : // if (movedist > tol * 100.0) nout=1;
1408 : : }
1409 : : }
1410 : :
1411 : : // Check if we are continuing to project outside the facet.
1412 : : // If so, then stop now
1413 : :
1414 [ + + ]: 8927 : if( nout > 3 ) { done = true; }
1415 : :
1416 : : // set up for next iteration
1417 : :
1418 [ + + ]: 8927 : if( !done )
1419 : : {
1420 [ + + ]: 8790 : if( distnew < bestdist )
1421 : : {
1422 : 5980 : bestdist = distnew;
1423 : 5980 : bestac = newac;
1424 : 5980 : bestpt = newpt;
1425 : 5980 : bestnorm = norm;
1426 : : }
1427 : 8790 : lastdist = distnew;
1428 : 8790 : lastpt = newpt;
1429 : 8790 : lastac = newac;
1430 [ + - ]: 11470 : move = pt - lastpt;
1431 : : }
1432 : : }
1433 : : }
1434 : :
1435 : 2680 : eval_pt = bestpt;
1436 [ - + ]: 2680 : if( eval_norm ) { *eval_norm = bestnorm; }
1437 : 2680 : outside = ( nout > 0 ) ? true : false;
1438 : 2680 : ac = bestac;
1439 : :
1440 : 2680 : return status;
1441 : : }
1442 : :
1443 : : //===========================================================================
1444 : : // Function Name: ac_at_edge
1445 : : //
1446 : : // Member Type: PRIVATE
1447 : : // Description: determine the area coordinate of the facet at the edge
1448 : : //===========================================================================
1449 : 0 : void SmoothFace::ac_at_edge( CartVect& fac, // facet area coordinate
1450 : : CartVect& eac, // edge area coordinate
1451 : : int edge_id ) // id of edge
1452 : : {
1453 : : double u, v, w;
1454 [ # # # # ]: 0 : switch( edge_id )
1455 : : {
1456 : : case 0:
1457 : 0 : u = 0.0;
1458 : 0 : v = fac[1] / ( fac[1] + fac[2] ); // v = fac.y() / (fac.y() + fac.z());
1459 : 0 : w = 1.0 - v;
1460 : 0 : break;
1461 : : case 1:
1462 : 0 : u = fac[0] / ( fac[0] + fac[2] ); // u = fac.x() / (fac.x() + fac.z());
1463 : 0 : v = 0.0;
1464 : 0 : w = 1.0 - u;
1465 : 0 : break;
1466 : : case 2:
1467 : 0 : u = fac[0] / ( fac[0] + fac[1] ); // u = fac.x() / (fac.x() + fac.y());
1468 : 0 : v = 1.0 - u;
1469 : 0 : w = 0.0;
1470 : 0 : break;
1471 : : default:
1472 : 0 : assert( 0 );
1473 : : u = -1; // needed to eliminate warnings about used before set
1474 : : v = -1; // needed to eliminate warnings about used before set
1475 : : w = -1; // needed to eliminate warnings about used before set
1476 : : break;
1477 : : }
1478 : 0 : eac[0] = u;
1479 : 0 : eac[1] = v;
1480 : 0 : eac[2] = w; //= CartVect(u, v, w);
1481 : 0 : }
1482 : :
1483 : : //===========================================================================
1484 : : // Function Name: project_to_facet
1485 : : //
1486 : : // Member Type: PUBLIC
1487 : : // Description: project to a single facet. Uses the input areacoord as
1488 : : // a starting guess.
1489 : : //===========================================================================
1490 : 2680 : ErrorCode SmoothFace::project_to_facet( EntityHandle facet, CartVect& pt, CartVect& areacoord, CartVect& close_point,
1491 : : bool& outside_facet, double compare_tol )
1492 : : {
1493 : 2680 : const EntityHandle* conn3 = NULL;
1494 : 2680 : int nnodes = 0;
1495 [ + - ]: 2680 : _mb->get_connectivity( facet, conn3, nnodes );
1496 : : //
1497 : : // double coords[9]; // store the coordinates for the nodes
1498 : : //_mb->get_coords(conn3, 3, coords);
1499 [ + - ][ + + ]: 10720 : CartVect p[3];
1500 [ + - ]: 2680 : _mb->get_coords( conn3, 3, (double*)&p[0] );
1501 : :
1502 : 2680 : int edge_id = -1;
1503 [ + - ]: 2680 : ErrorCode stat = project_to_patch( facet, areacoord, pt, close_point, NULL, outside_facet, compare_tol, edge_id );
1504 : : /* }
1505 : : break;
1506 : : }*/
1507 : :
1508 : 2680 : return stat;
1509 : : }
1510 : :
1511 : : //===========================================================================
1512 : : // Function Name: is_at_vertex
1513 : : //
1514 : : // Member Type: PRIVATE
1515 : : // Description: determine if the point is at one of the facet's vertices
1516 : : //===========================================================================
1517 : 2680 : bool SmoothFace::is_at_vertex( EntityHandle facet, // (IN) facet we are evaluating
1518 : : CartVect& pt, // (IN) the point
1519 : : CartVect& ac, // (IN) the ac of the point on the facet plane
1520 : : double compare_tol, // (IN) return TRUE of closer than this
1521 : : CartVect& eval_pt, // (OUT) location at vertex if TRUE
1522 : : CartVect* eval_norm_ptr ) // (OUT) normal at vertex if TRUE
1523 : : {
1524 : : double dist;
1525 [ + - ]: 2680 : CartVect vert_loc;
1526 : 2680 : const double actol = 0.1;
1527 : : // get coordinates get_coords
1528 : 2680 : const EntityHandle* conn3 = NULL;
1529 : 2680 : int nnodes = 0;
1530 [ + - ]: 2680 : _mb->get_connectivity( facet, conn3, nnodes );
1531 : : //
1532 : : // double coords[9]; // store the coordinates for the nodes
1533 : : //_mb->get_coords(conn3, 3, coords);
1534 [ + - ][ + + ]: 10720 : CartVect p[3];
1535 [ + - ]: 2680 : _mb->get_coords( conn3, 3, (double*)&p[0] );
1536 : : // also get the normals at nodes
1537 [ + - ][ + + ]: 10720 : CartVect NN[3];
1538 [ + - ]: 2680 : _mb->tag_get_data( _gradientTag, conn3, 3, (double*)&NN[0] );
1539 [ + - ][ + + ]: 2680 : if( fabs( ac[0] ) < actol && fabs( ac[1] ) < actol )
[ + - ][ + + ]
[ + + ]
1540 : : {
1541 : 154 : vert_loc = p[2];
1542 [ + - ][ + - ]: 154 : dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
1543 [ - + ]: 154 : if( dist <= compare_tol )
1544 : : {
1545 : 0 : eval_pt = vert_loc;
1546 [ # # ]: 0 : if( eval_norm_ptr ) { *eval_norm_ptr = NN[2]; }
1547 : 0 : return true;
1548 : : }
1549 : : }
1550 : :
1551 [ + - ][ + + ]: 2680 : if( fabs( ac[0] ) < actol && fabs( ac[2] ) < actol )
[ + - ][ + + ]
[ + + ]
1552 : : {
1553 : 145 : vert_loc = p[1];
1554 [ + - ][ + - ]: 145 : dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
1555 [ - + ]: 145 : if( dist <= compare_tol )
1556 : : {
1557 : 0 : eval_pt = vert_loc;
1558 [ # # ]: 0 : if( eval_norm_ptr )
1559 : : {
1560 : 0 : *eval_norm_ptr = NN[1]; // facet->point(1)->normal( facet );
1561 : : }
1562 : 0 : return true;
1563 : : }
1564 : : }
1565 : :
1566 [ + - ][ + + ]: 2680 : if( fabs( ac[1] ) < actol && fabs( ac[2] ) < actol )
[ + - ][ + + ]
[ + + ]
1567 : : {
1568 : 129 : vert_loc = p[0];
1569 [ + - ][ + - ]: 129 : dist = ( pt - vert_loc ).length(); // pt.distance_between( vert_loc );
1570 [ - + ]: 129 : if( dist <= compare_tol )
1571 : : {
1572 : 0 : eval_pt = vert_loc;
1573 [ # # ]: 0 : if( eval_norm_ptr ) { *eval_norm_ptr = NN[0]; }
1574 : 0 : return true;
1575 : : }
1576 : : }
1577 : :
1578 : 2680 : return false;
1579 : : }
1580 : :
1581 : : //===========================================================================
1582 : : // Function Name: move_ac_inside
1583 : : //
1584 : : // Member Type: PRIVATE
1585 : : // Description: find the closest area coordinate to the boundary of the
1586 : : // patch if any of its components are < 0
1587 : : // Return if the ac was modified.
1588 : : //===========================================================================
1589 : 4790 : bool SmoothFace::move_ac_inside( CartVect& ac, double tol )
1590 : : {
1591 : 4790 : int nout = 0;
1592 [ + + ]: 4790 : if( ac[0] < -tol )
1593 : : {
1594 : 600 : ac[0] = 0.0;
1595 : 600 : ac[1] = ac[1] / ( ac[1] + ac[2] ); //( ac.y() / (ac.y() + ac.z()) ;
1596 : 600 : ac[2] = 1. - ac[1]; // ac.z( 1.0 - ac.y() );
1597 : 600 : nout++;
1598 : : }
1599 [ + + ]: 4790 : if( ac[1] < -tol )
1600 : : {
1601 : 794 : ac[1] = 0.; // ac.y( 0.0 );
1602 : 794 : ac[0] = ac[0] / ( ac[0] + ac[2] ); // ac.x( ac.x() / (ac.x() + ac.z()) );
1603 : 794 : ac[2] = 1. - ac[0]; // ac.z( 1.0 - ac.x() );
1604 : 794 : nout++;
1605 : : }
1606 [ + + ]: 4790 : if( ac[2] < -tol )
1607 : : {
1608 : 732 : ac[2] = 0.; // ac.z( 0.0 );
1609 : 732 : ac[0] = ac[0] / ( ac[0] + ac[1] ); // ac.x( ac.x() / (ac.x() + ac.y()) );
1610 : 732 : ac[1] = 1. - ac[0]; // ac.y( 1.0 - ac.x() );
1611 : 732 : nout++;
1612 : : }
1613 : 4790 : return ( nout > 0 ) ? true : false;
1614 : : }
1615 : :
1616 : : //===========================================================================
1617 : : // Function Name: hodograph
1618 : : //
1619 : : // Member Type: PUBLIC
1620 : : // Description: get the hodograph control points for the facet
1621 : : // Note: This is a triangle cubic patch that is defined by the
1622 : : // normals of quartic facet control point lattice. Returned coordinates
1623 : : // in Nijk are defined by the following diagram
1624 : : //
1625 : : //
1626 : : // *9 index polar
1627 : : // / \ 0 300 point(0)
1628 : : // / \ 1 210
1629 : : // 7*-----*8 2 120
1630 : : // / \ / \ 3 030 point(1)
1631 : : // / \ / \ 4 201
1632 : : // 4*----5*-----*6 5 111
1633 : : // / \ / \ / \ 6 021
1634 : : // / \ / \ / \ 7 102
1635 : : // *-----*-----*-----* 8 012
1636 : : // 0 1 2 3 9 003 point(2)
1637 : : //
1638 : :
1639 : : //===========================================================================
1640 : : // Function Name: eval_bezier_patch_normal
1641 : : //
1642 : : // Member Type: PRIVATE
1643 : : // Description: evaluate the Bezier patch defined at a facet
1644 : : //===========================================================================
1645 : 14 : ErrorCode SmoothFace::eval_bezier_patch_normal( EntityHandle facet, CartVect& areacoord, CartVect& normal )
1646 : : {
1647 : : // interpolate internal control points
1648 : :
1649 [ + - ][ + + ]: 98 : CartVect gctrl_pts[6];
1650 : : // facet->get_control_points( gctrl_pts );
1651 [ + - ]: 14 : ErrorCode rval = _mb->tag_get_data( _facetCtrlTag, &facet, 1, &gctrl_pts[0] );
1652 [ - + ]: 14 : assert( rval == MB_SUCCESS );
1653 [ - + ]: 14 : if( MB_SUCCESS != rval ) return rval;
1654 : : // _gradientTag
1655 : : // get normals at points
1656 : 14 : const EntityHandle* conn3 = NULL;
1657 : 14 : int nnodes = 0;
1658 [ + - ]: 14 : rval = _mb->get_connectivity( facet, conn3, nnodes );
1659 [ - + ]: 14 : if( MB_SUCCESS != rval ) return rval;
1660 : :
1661 [ + - ][ + + ]: 56 : CartVect NN[3];
1662 [ + - ]: 14 : rval = _mb->tag_get_data( _gradientTag, conn3, 3, &NN[0] );
1663 [ - + ]: 14 : assert( rval == MB_SUCCESS );
1664 [ - + ]: 14 : if( MB_SUCCESS != rval ) return rval;
1665 : :
1666 [ + - ][ + - ]: 14 : if( fabs( areacoord[1] + areacoord[2] ) < 1.0e-6 )
[ - + ]
1667 : : {
1668 : 0 : normal = NN[0];
1669 : 0 : return MB_SUCCESS;
1670 : : }
1671 [ + - ][ + - ]: 14 : if( fabs( areacoord[0] + areacoord[2] ) < 1.0e-6 )
[ - + ]
1672 : : {
1673 : 0 : normal = NN[1]; // facet->point(1)->normal(facet);
1674 : 0 : return MB_SUCCESS;
1675 : : }
1676 [ + - ][ + - ]: 14 : if( fabs( areacoord[0] + areacoord[1] ) < 1.0e-6 )
[ + + ]
1677 : : {
1678 : 1 : normal = NN[2]; // facet->point(2)->normal(facet);
1679 : 1 : return MB_SUCCESS;
1680 : : }
1681 : :
1682 : : // compute the hodograph of the quartic Gregory patch
1683 : :
1684 [ + - ][ + + ]: 143 : CartVect Nijk[10];
1685 : : // hodograph(facet,areacoord,Nijk);
1686 : : // start copy from hodograph
1687 : : // CubitVector gctrl_pts[6];
1688 : : // facet->get_control_points( gctrl_pts );
1689 [ + - ][ + + ]: 52 : CartVect P_facet[3];
1690 : :
1691 : : // 2,1,1
1692 : : /*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
1693 : : (areacoord.y() * gctrl_pts[3] +
1694 : : areacoord.z() * gctrl_pts[4]);*/
1695 : : P_facet[0] =
1696 [ + - ][ + - ]: 13 : ( 1.0e0 / ( areacoord[1] + areacoord[2] ) ) * ( areacoord[1] * gctrl_pts[3] + areacoord[2] * gctrl_pts[4] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1697 : : // 1,2,1
1698 : : /*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
1699 : : (areacoord.x() * gctrl_pts[0] +
1700 : : areacoord.z() * gctrl_pts[5]);*/
1701 : : P_facet[1] =
1702 [ + - ][ + - ]: 13 : ( 1.0e0 / ( areacoord[0] + areacoord[2] ) ) * ( areacoord[0] * gctrl_pts[0] + areacoord[2] * gctrl_pts[5] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1703 : : // 1,1,2
1704 : : /*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
1705 : : (areacoord.x() * gctrl_pts[1] +
1706 : : areacoord.y() * gctrl_pts[2]);*/
1707 : : P_facet[2] =
1708 [ + - ][ + - ]: 13 : ( 1.0e0 / ( areacoord[0] + areacoord[1] ) ) * ( areacoord[0] * gctrl_pts[1] + areacoord[1] * gctrl_pts[2] );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1709 : :
1710 : : // corner control points are just the normals at the points
1711 : :
1712 : : // 3, 0, 0
1713 : 13 : Nijk[0] = NN[0];
1714 : : // 0, 3, 0
1715 : 13 : Nijk[3] = NN[1];
1716 : : // 0, 0, 3
1717 : 13 : Nijk[9] = NN[2]; // facet->point(2)->normal(facet);
1718 : :
1719 : : // fill in the boundary control points. Define as the normal to the local
1720 : : // triangle formed by the quartic control point lattice
1721 : :
1722 : : // store here again the 9 control points on the edges
1723 [ + - ][ + + ]: 130 : CartVect CP[9]; // 9 control points on the edges,
1724 [ + - ]: 13 : rval = _mb->tag_get_data( _facetEdgeCtrlTag, &facet, 1, &CP[0] );
1725 [ - + ]: 13 : if( MB_SUCCESS != rval ) return rval;
1726 : : // there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
1727 : : // CubitFacetEdge *edge;
1728 : : // edge = facet->edge( 2 );
1729 : : // CubitVector ctrl_pts[5];
1730 : : // edge->control_points(facet, ctrl_pts);
1731 : :
1732 : : // 2, 1, 0
1733 : : // Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
1734 [ + - ][ + - ]: 13 : Nijk[1] = ( CP[7] - CP[6] ) * ( P_facet[0] - CP[6] );
[ + - ]
1735 [ + - ]: 13 : Nijk[1].normalize();
1736 : :
1737 : : // 1, 2, 0
1738 : : // Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
1739 [ + - ][ + - ]: 13 : Nijk[2] = ( CP[8] - CP[7] ) * ( P_facet[1] - CP[7] );
[ + - ]
1740 [ + - ]: 13 : Nijk[2].normalize();
1741 : :
1742 : : // edge = facet->edge( 0 );
1743 : : // edge->control_points(facet, ctrl_pts);
1744 : :
1745 : : // 0, 2, 1
1746 : : // Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
1747 [ + - ][ + - ]: 13 : Nijk[6] = ( CP[0] - P_facet[1] ) * ( CP[1] - P_facet[1] );
[ + - ]
1748 [ + - ]: 13 : Nijk[6].normalize();
1749 : :
1750 : : // 0, 1, 2
1751 : : // Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
1752 [ + - ][ + - ]: 13 : Nijk[8] = ( CP[1] - P_facet[2] ) * ( CP[2] - P_facet[2] );
[ + - ]
1753 [ + - ]: 13 : Nijk[8].normalize();
1754 : :
1755 : : // edge = facet->edge( 1 );
1756 : : // edge->control_points(facet, ctrl_pts);
1757 : :
1758 : : // 1, 0, 2
1759 : : // Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
1760 [ + - ][ + - ]: 13 : Nijk[7] = ( P_facet[2] - CP[4] ) * ( CP[3] - CP[4] );
[ + - ]
1761 [ + - ]: 13 : Nijk[7].normalize();
1762 : :
1763 : : // 2, 0, 1
1764 : : // Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
1765 [ + - ][ + - ]: 13 : Nijk[4] = ( P_facet[0] - CP[5] ) * ( CP[4] - CP[5] );
[ + - ]
1766 [ + - ]: 13 : Nijk[4].normalize();
1767 : :
1768 : : // 1, 1, 1
1769 [ + - ][ + - ]: 13 : Nijk[5] = ( P_facet[1] - P_facet[0] ) * ( P_facet[2] - P_facet[0] );
[ + - ]
1770 [ + - ]: 13 : Nijk[5].normalize();
1771 : : // end copy from hodograph
1772 : :
1773 : : // sum the contribution from each of the control points
1774 : :
1775 [ + - ]: 13 : normal = CartVect( 0.0e0, 0.0e0, 0.0e0 );
1776 : :
1777 : : // i=3; j=0; k=0;
1778 : : // double Bsum = 0.0;
1779 [ + - ][ + - ]: 13 : double B = mbcube( areacoord[0] );
[ + - ]
1780 : : // Bsum += B;
1781 [ + - ][ + - ]: 13 : normal += B * Nijk[0];
1782 : :
1783 : : // i=2; j=1; k=0;
1784 [ + - ][ + - ]: 13 : B = 3.0 * mbsqr( areacoord[0] ) * areacoord[1];
[ + - ]
1785 : : // Bsum += B;
1786 [ + - ][ + - ]: 13 : normal += B * Nijk[1];
1787 : :
1788 : : // i=1; j=2; k=0;
1789 [ + - ][ + - ]: 13 : B = 3.0 * areacoord[0] * mbsqr( areacoord[1] );
[ + - ]
1790 : : // Bsum += B;
1791 [ + - ][ + - ]: 13 : normal += B * Nijk[2];
1792 : :
1793 : : // i=0; j=3; k=0;
1794 [ + - ][ + - ]: 13 : B = mbcube( areacoord[1] );
[ + - ]
1795 : : // Bsum += B;
1796 [ + - ][ + - ]: 13 : normal += B * Nijk[3];
1797 : :
1798 : : // i=2; j=0; k=1;
1799 [ + - ][ + - ]: 13 : B = 3.0 * mbsqr( areacoord[0] ) * areacoord[2];
[ + - ]
1800 : : // Bsum += B;
1801 [ + - ][ + - ]: 13 : normal += B * Nijk[4];
1802 : :
1803 : : // i=1; j=1; k=1;
1804 [ + - ][ + - ]: 13 : B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
[ + - ]
1805 : : // Bsum += B;
1806 [ + - ][ + - ]: 13 : normal += B * Nijk[5];
1807 : :
1808 : : // i=0; j=2; k=1;
1809 [ + - ][ + - ]: 13 : B = 3.0 * mbsqr( areacoord[1] ) * areacoord[2];
[ + - ]
1810 : : // Bsum += B;
1811 [ + - ][ + - ]: 13 : normal += B * Nijk[6];
1812 : :
1813 : : // i=1; j=0; k=2;
1814 [ + - ][ + - ]: 13 : B = 3.0 * areacoord[0] * mbsqr( areacoord[2] );
[ + - ]
1815 : : // Bsum += B;
1816 [ + - ][ + - ]: 13 : normal += B * Nijk[7];
1817 : :
1818 : : // i=0; j=1; k=2;
1819 [ + - ][ + - ]: 13 : B = 3.0 * areacoord[1] * mbsqr( areacoord[2] );
[ + - ]
1820 : : // Bsum += B;
1821 [ + - ][ + - ]: 13 : normal += B * Nijk[8];
1822 : :
1823 : : // i=0; j=0; k=3;
1824 [ + - ][ + - ]: 13 : B = mbcube( areacoord[2] );
[ + - ]
1825 : : // Bsum += B;
1826 [ + - ][ + - ]: 13 : normal += B * Nijk[9];
1827 : :
1828 : : // assert(fabs(Bsum - 1.0) < 1e-9);
1829 : :
1830 [ + - ]: 13 : normal.normalize();
1831 : :
1832 : 14 : return MB_SUCCESS;
1833 : : }
1834 : :
1835 : 2667 : ErrorCode SmoothFace::get_normals_for_vertices( const EntityHandle* conn2, CartVect N[2] )
1836 : : // this method will be called to retrieve the normals needed in the calculation of control edge
1837 : : // points..
1838 : : {
1839 : : // CartVect N[2];
1840 : : //_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
1841 : 2667 : ErrorCode rval = _mb->tag_get_data( _gradientTag, conn2, 2, (double*)&N[0] );
1842 : 2667 : return rval;
1843 : : }
1844 : :
1845 : 4 : ErrorCode SmoothFace::ray_intersection_correct( EntityHandle, // (IN) the facet where the patch is defined
1846 : : CartVect& pt, // (IN) shoot from
1847 : : CartVect& ray, // (IN) ray direction
1848 : : CartVect& eval_pt, // (INOUT) The intersection point
1849 : : double& distance, // (IN OUT) the new distance
1850 : : bool& outside )
1851 : : {
1852 : : // find a point on the smooth surface
1853 : 4 : CartVect currentPoint = eval_pt;
1854 : 4 : int numIter = 0;
1855 : 4 : double improvement = 1.e20;
1856 [ + - ]: 4 : CartVect diff;
1857 [ + - ][ + + ]: 10 : while( numIter++ < 5 && improvement > 0.01 )
[ + + ]
1858 : : {
1859 [ + - ]: 6 : CartVect newPos;
1860 : :
1861 : 6 : bool trim = false; // is it needed?
1862 : 6 : outside = true;
1863 [ + - ]: 6 : CartVect closestPoint;
1864 [ + - ]: 6 : CartVect normal;
1865 : :
1866 [ + - ]: 6 : ErrorCode rval = project_to_facets_main( currentPoint, trim, outside, &newPos, &normal );
1867 [ - + ]: 6 : if( MB_SUCCESS != rval ) return rval;
1868 [ - + ]: 6 : assert( rval == MB_SUCCESS );
1869 [ + - ]: 6 : diff = newPos - currentPoint;
1870 [ + - ]: 6 : improvement = diff.length();
1871 : : // ( pt + t * ray - closest ) % normal = 0;
1872 : : // intersect tangent plane that goes through closest point with the direction
1873 : : // t = normal%(closest-pt) / normal%ray;
1874 [ + - ]: 6 : double dot = normal % ray; // if it is 0, get out while we can
1875 [ - + ]: 6 : if( dot < 0.00001 )
1876 : : {
1877 : : // bad convergence, get out, do not modify anything
1878 : 0 : return MB_SUCCESS;
1879 : : }
1880 [ + - ][ + - ]: 6 : double t = ( ( newPos - pt ) % normal ) / ( dot );
1881 [ + - ][ + - ]: 6 : currentPoint = pt + t * ray;
1882 : : }
1883 : 4 : eval_pt = currentPoint;
1884 [ + - ]: 4 : diff = currentPoint - pt;
1885 [ + - ]: 4 : distance = diff.length();
1886 : 4 : return MB_SUCCESS;
1887 : : }
1888 [ + - ][ + - ]: 44 : } // namespace moab
|
