Branch data Line data Source code
1 : : //- Class: FacetEvalTool
2 : : //- Description: The FacetEvalTool is a tool to perform generic geometric
3 : : //- operations on a set of facets.
4 : : //- Assumptions: All of the facets are connected topologically correct.
5 : : //-
6 : : //- Owner: Steve J. Owen
7 : : //- Checked by:
8 : :
9 : : #include <float.h>
10 : : #include "CastTo.hpp"
11 : : #include "CubitVector.hpp"
12 : : #include "CubitPoint.hpp"
13 : : #include "CubitFacet.hpp"
14 : : #include "CubitFacetEdge.hpp"
15 : : #include "CubitFacetEdgeData.hpp"
16 : : #include "CpuTimer.hpp"
17 : : #include "CubitMessage.hpp"
18 : : #include "CubitBox.hpp"
19 : : #include "DLIList.hpp"
20 : : #include "FacetEvalTool.hpp"
21 : : #include "GeometryDefines.h"
22 : : #include "CubitTransformMatrix.hpp"
23 : : #include "ChollaDebug.hpp"
24 : : #include "TDFacetBoundaryEdge.hpp"
25 : : #include "TDFacetBoundaryPoint.hpp"
26 : : #include "GfxDebug.hpp"
27 : : //#include "FacetParamTool.hpp"
28 : : #include "KDDTree.hpp"
29 : : #include "AbstractTree.hpp"
30 : : #include "CubitFileIOWrapper.hpp"
31 : : #include "FacetDataUtil.hpp"
32 : :
33 : : double FacetEvalTool::timeGridSearch = 0.0;
34 : : double FacetEvalTool::timeFacetProject = 0.0;
35 : : int FacetEvalTool::numEvals = 0;
36 : : #define GRID_SEARCH_THRESHOLD 20
37 : :
38 : : //===========================================================================
39 : : //Function Name: FacetEvalTool
40 : : //
41 : : //Member Type: PUBLIC
42 : : //Description: constructor
43 : : //===========================================================================
44 : 341 : CubitStatus FacetEvalTool::initialize(
45 : : DLIList<CubitFacet*> &facet_list,
46 : : DLIList<CubitPoint*> &point_list, // if this isn't filled in, I'll do it here
47 : : int interp_order, // 0 = linear or 4 = Bezier only for now
48 : : double min_dot ) // from feature angle
49 : : {
50 : 341 : myFacetList = facet_list;
51 : 341 : myPointList = point_list;
52 : 341 : isParameterized = CUBIT_FALSE;
53 : 341 : minDot = min_dot;
54 : 341 : output_id = -1;
55 : :
56 : : // mark all of the facets on this surface with the toolID
57 : :
58 : : int i;
59 [ + + ]: 1397 : for (i = 0; i< myFacetList.size(); i++)
60 : 1056 : myFacetList.get_and_step()->set_tool_id( toolID );
61 : :
62 : : // generate the list of points, if not already defined
63 : :
64 [ - + ]: 341 : if (myPointList.size() == 0)
65 : : {
66 : 0 : get_points_from_facets( facet_list, myPointList );
67 : : }
68 : 341 : interpOrder = 0; // default to project to linear facet
69 : :
70 : : // set the bounding box and compareTol and setup grid search
71 : 341 : myBBox = NULL;
72 : 341 : aTree = NULL;
73 : :
74 : 341 : reset_bounding_box();
75 : :
76 : :
77 : 341 : lastFacet = NULL;
78 : 341 : isFlat = -999;
79 : 341 : isFlat = is_flat();
80 : :
81 [ - + ]: 341 : if (interp_order == -1) {
82 : 0 : interpOrder = 0; // use linear as default interp order
83 : : }
84 : : else {
85 : 341 : interpOrder = interp_order;
86 : : }
87 : :
88 : : // generate edges
89 : :
90 [ + - - + ]: 682 : if ( CUBIT_SUCCESS != get_edges_from_facets( facet_list, myEdgeList ) ||
[ - + ]
91 : 341 : CUBIT_SUCCESS != get_loops_from_facets( myEdgeList, myLoopList ) )
92 : : {
93 [ # # ][ # # ]: 0 : PRINT_ERROR("Unable to initialize Facet Evaluation Tool.\n");
94 : 0 : return CUBIT_FAILURE;
95 : : }
96 : :
97 [ + + ]: 341 : if (interpOrder > 0) {
98 : 66 : init_gradient();
99 [ - + ]: 66 : if (interpOrder == 3) { // least_squares
100 : 0 : init_quadrics();
101 : : }
102 [ + - ]: 66 : else if(interpOrder == 4) { // spline
103 : 66 : init_bezier_surface();
104 : : }
105 : : }
106 : :
107 : : // compute the area of all facets
108 : :
109 : 341 : myArea = -1.0;
110 : 341 : myArea = area();
111 : :
112 : :
113 : 341 : int mydebug = 0;
114 [ - + ]: 341 : if (mydebug)
115 : : {
116 : 0 : debug_draw_eval_bezier_facet( myFacetList.get_and_step() );
117 : 0 : debug_draw_facets();
118 : 0 : debug_draw_facet_normals();
119 [ # # ]: 0 : if (interpOrder > 0) debug_draw_point_normals();
120 : 0 : GfxDebug::flush();
121 : 0 : GfxDebug::mouse_xforms();
122 : : }
123 : :
124 : : //parameterize();
125 : 341 : return CUBIT_SUCCESS;
126 : : }
127 : :
128 : : //===========================================================================
129 : : //Function Name: FacetEvalTool
130 : : //Member Type: PUBLIC
131 : : //Descriptoin: constructor. This one is an empty constructor
132 : : // sed with restore method
133 : : //===========================================================================
134 [ + - ][ + - ]: 946 : FacetEvalTool::FacetEvalTool()
[ + - ][ + - ]
135 : : {
136 : : static int counter = 1;
137 : 473 : toolID = counter++;
138 : 473 : myBBox = NULL;
139 : 473 : aTree = NULL;
140 : 473 : lastFacet = NULL;
141 : 473 : isFlat = -999;
142 : 473 : myArea = -1.0;
143 : 473 : interpOrder = 0;
144 : 473 : minDot = 0;
145 : 473 : isParameterized = CUBIT_FALSE;
146 : 473 : }
147 : :
148 : : //===========================================================================
149 : : //Function Name: ~FacetEvalTool
150 : : //
151 : : //Member Type: PUBLIC
152 : : //Descriptoin: destructor
153 : : //===========================================================================
154 [ # # ][ # # ]: 0 : FacetEvalTool::~FacetEvalTool()
[ # # ][ # # ]
155 : : {
156 [ # # ][ # # ]: 0 : if ( DEBUG_FLAG(110) )
[ # # ]
157 : : {
158 [ # # ][ # # ]: 0 : PRINT_INFO("num facets = %d\n",myFacetList.size());
[ # # ][ # # ]
[ # # ]
159 [ # # ][ # # ]: 0 : PRINT_INFO("numEvals = %d\n",numEvals);
[ # # ][ # # ]
160 [ # # ][ # # ]: 0 : PRINT_INFO("timeFacetProject = %f\n",timeFacetProject);
[ # # ][ # # ]
161 [ # # ][ # # ]: 0 : PRINT_INFO("timeGridSearch = %f\n",timeGridSearch);
[ # # ][ # # ]
162 : : }
163 : :
164 [ # # ]: 0 : if (aTree != NULL)
165 [ # # ][ # # ]: 0 : delete aTree;
166 : :
167 [ # # ][ # # ]: 0 : if (myBBox) delete myBBox;
[ # # ]
168 [ # # ]: 0 : destroy_facets();
169 : :
170 [ # # ]: 0 : myLoopList.reset();
171 : : int i;
172 [ # # ][ # # ]: 0 : for (i=myLoopList.size(); i>0; i--)
173 : : {
174 [ # # ][ # # ]: 0 : delete myLoopList.get_and_step();
[ # # ]
175 : : }
176 [ # # ]: 0 : myLoopList.clean_out();
177 : 0 : }
178 : :
179 : 0 : void FacetEvalTool::remove_facets( DLIList<CubitFacet*>& facets )
180 : : {
181 : 0 : facets = myFacetList;
182 [ # # ]: 0 : for ( int i = facets.size(); i--; )
183 : 0 : facets.step_and_get()->set_tool_id(0);
184 : 0 : myFacetList.clean_out();
185 : 0 : myEdgeList.clean_out();
186 : 0 : myPointList.clean_out();
187 : 0 : }
188 : :
189 : 0 : CubitStatus FacetEvalTool::reverse_facets( )
190 : : {
191 : : int i;
192 : : CubitFacet* temp_facet;
193 : :
194 [ # # ]: 0 : for(i=0; i<myFacetList.size(); i++){
195 : 0 : temp_facet=myFacetList.get_and_step();
196 [ # # ]: 0 : if(!temp_facet){
197 [ # # ][ # # ]: 0 : PRINT_ERROR("Unexpected NULL pointer for facet.\n");
198 : 0 : return CUBIT_FAILURE;
199 : : }
200 : 0 : temp_facet->flip();
201 : : }
202 : :
203 : 0 : return CUBIT_SUCCESS;
204 : : }
205 : :
206 : 0 : CubitStatus FacetEvalTool::reverse_facets( DLIList<CubitFacet *> &facets )
207 : : {
208 : : int i;
209 : : CubitFacet* temp_facet;
210 : :
211 [ # # ]: 0 : for(i=0; i<facets.size(); i++){
212 : 0 : temp_facet=facets.get_and_step();
213 [ # # ]: 0 : if(!temp_facet){
214 [ # # ][ # # ]: 0 : PRINT_ERROR("Unexpected NULL pointer for facet.\n");
215 : 0 : return CUBIT_FAILURE;
216 : : }
217 : 0 : temp_facet->flip();
218 : : }
219 : :
220 : 0 : return CUBIT_SUCCESS;
221 : : }
222 : :
223 : :
224 : 0 : void FacetEvalTool::replace_facets( DLIList<CubitFacet *> &facet_list )
225 : : {
226 : 0 : myFacetList.clean_out();
227 : 0 : myEdgeList.clean_out();
228 : 0 : myPointList.clean_out();
229 : 0 : FacetDataUtil::get_points( facet_list, myPointList );
230 : 0 : FacetDataUtil::get_edges( facet_list, myEdgeList );
231 : 0 : myFacetList = facet_list;
232 [ # # ]: 0 : for ( int i = myFacetList.size(); i--; )
233 : 0 : myFacetList.step_and_get()->set_tool_id(toolID);
234 : 0 : reset_bounding_box();
235 : 0 : }
236 : :
237 : :
238 : : //===========================================================================
239 : : //Function Name: set_up_grid_search
240 : : //
241 : : //Member Type: PUBLIC
242 : : //Descriptoin: set up grid search if we have a lot of facets
243 : : //===========================================================================
244 : 473 : void FacetEvalTool::set_up_grid_search(
245 : : double geom_factor )
246 : : {
247 [ - + ]: 473 : if(aTree)
248 [ # # ]: 0 : delete aTree;
249 : :
250 [ + - ]: 473 : if (myFacetList.size() < GRID_SEARCH_THRESHOLD)
251 : : {
252 : 473 : aTree = NULL;
253 : : }
254 : : else
255 : : {
256 [ # # ]: 0 : aTree = new KDDTree<CubitFacet*> (GEOMETRY_RESABS*geom_factor, false/*self-balancing off*/);
257 [ # # ]: 0 : for ( int ii = myFacetList.size(); ii > 0; ii-- )
258 : : {
259 : 0 : aTree->add(myFacetList.get_and_step());
260 : : }
261 : 0 : aTree->balance();
262 : : }
263 : 473 : }
264 : :
265 : : //===========================================================================
266 : : //Function Name: is_flat
267 : : //
268 : : //Member Type: PUBLIC
269 : : //Descriptoin: Determine if the facets are all in the same plane
270 : : //===========================================================================
271 : 1221 : int FacetEvalTool::is_flat()
272 : : {
273 : : int ii;
274 [ + + ]: 1221 : if (isFlat != -999) {
275 : 880 : return isFlat;
276 : : }
277 : : else {
278 : 341 : isFlat = CUBIT_TRUE;
279 [ + - ][ + - ]: 341 : CubitVector firstnormal = myFacetList.get_and_step()->normal();
280 [ + - ]: 341 : firstnormal.normalize();
281 [ + - ]: 341 : CubitVector normal;
282 [ + - ][ + + ]: 1034 : for ( ii = myFacetList.size(); ii > 1 &&isFlat; ii-- ){
[ + + ]
283 [ + - ][ + - ]: 693 : normal = myFacetList.get_and_step()->normal();
[ + - ]
284 [ + - ]: 693 : normal.normalize();
285 [ + - ][ + + ]: 693 : if (fabs(normal%firstnormal) < 0.9987654321) {
286 : 33 : isFlat = CUBIT_FALSE;
287 : : }
288 : : }
289 : : }
290 : 1221 : return isFlat;
291 : : }
292 : :
293 : : //===========================================================================
294 : : //Function Name: init_gradient
295 : : //
296 : : //Member Type: PRIVATE
297 : : //Descriptoin: initialize the gradients for order 1 interpolation
298 : : //===========================================================================
299 : 66 : CubitStatus FacetEvalTool::init_gradient()
300 : : {
301 : : int i,j;
302 : :
303 : : // retrieve all faces attached to the points in point_list
304 : :
305 [ + + ]: 330 : for (i = 0; i < myPointList.size(); i++)
306 : : {
307 [ + - ]: 264 : CubitPoint* point = myPointList.get_and_step();
308 : :
309 [ + - ]: 264 : DLIList<CubitFacet*> adj_facet_list;
310 [ + - ]: 264 : point->facets(adj_facet_list);
311 [ + - ][ + - ]: 264 : if (adj_facet_list.size() > 0) {
312 [ + - ]: 264 : CubitVector avg_normal(0.0e0, 0.0e0, 0.0e0);
313 : 264 : double totangle = 0.0e0;
314 : :
315 : : // weight the normal by the spanning angle at the point
316 : :
317 [ + - ][ + + ]: 1452 : for (j = 0; j < adj_facet_list.size(); j++)
318 : : {
319 [ + - ]: 1188 : CubitFacet* facet = adj_facet_list.get_and_step();
320 [ + - ]: 1188 : double angle = facet->angle( point );
321 [ + - ]: 1188 : facet->weight( angle );
322 : 1188 : totangle += angle;
323 : : }
324 [ + - ][ + + ]: 1452 : for (j = 0; j < adj_facet_list.size(); j++)
325 : : {
326 [ + - ]: 1188 : CubitFacet* facet = adj_facet_list.get_and_step();
327 [ + - ]: 1188 : CubitVector normal = facet->normal();
328 [ + - ]: 1188 : normal.normalize();
329 [ + - ][ + - ]: 1188 : avg_normal += (facet->weight() / totangle) * normal;
[ + - ]
330 : : }
331 [ + - ]: 264 : avg_normal.normalize();
332 [ + - ]: 264 : point->normal(avg_normal);
333 [ + - ][ + - ]: 264 : double coefd = -(point->coordinates()%avg_normal);
334 [ + - ]: 264 : point->d_coef( coefd );
335 : : }
336 [ + - ]: 264 : }
337 : 66 : return CUBIT_SUCCESS;
338 : : }
339 : :
340 : : //===========================================================================
341 : : //Function Name: init_quadrics
342 : : // NOTE: I (Roshan) fixed couple of bugs on Aug 02, 2010. See BUGFIX comments below. I didn't test this function; however, I have tested CMLSmoothTool::init_quadric().
343 : : //Member Type: PRIVATE
344 : : //Descriptoin: initialize the quadrics at the facet vertices for order 2
345 : : // interpolation
346 : : //===========================================================================
347 : 0 : CubitStatus FacetEvalTool::init_quadrics()
348 : : {
349 : : // use the normal at the vertices as a local space with which to do
350 : : // interpolation. A quadric will be approximated from the surrounding
351 : : // vertices. Interpolation will provide a "z" deviation from the tangent
352 : : // plane at the vertex
353 : :
354 : : // define a basis set of vectors at each point (assumes the gradients
355 : : // have already been approximated
356 : :
357 : : int i,j;
358 : : CubitPoint* point;
359 [ # # ][ # # ]: 0 : for (i = 0; i < myPointList.size(); i++)
360 : : {
361 [ # # ]: 0 : point = myPointList.get_and_step();
362 [ # # ]: 0 : point->define_tangent_vectors();
363 : : }
364 : :
365 : : // set up for least squares
366 : :
367 : : CubitStatus status;
368 : : #define MAX_CLOSE_POINTS 100
369 : : CubitPoint *close_points[MAX_CLOSE_POINTS];
370 [ # # ][ # # ]: 0 : CubitVector coords[MAX_CLOSE_POINTS], cp;
[ # # ]
371 : : double weight[MAX_CLOSE_POINTS];
372 : : int num_close;
373 [ # # ]: 0 : CubitMatrix lhs(5,5);
374 [ # # ][ # # ]: 0 : CubitMatrix rhs(5,1);
375 [ # # ][ # # ]: 0 : CubitMatrix coef(5,1);
376 [ # # ][ # # ]: 0 : for (i = 0; i < myPointList.size(); i++)
377 : : {
378 [ # # ]: 0 : point = myPointList.get_and_step();
379 : : status = get_close_points( point, close_points, num_close,
380 [ # # ]: 0 : MAX_CLOSE_POINTS, 5 );
381 [ # # ]: 0 : if (status != CUBIT_SUCCESS) {
382 : 0 : return status;
383 : : }
384 : :
385 : : // transform to local system in x-y
386 : : // determine weights based on inverse distance
387 : :
388 : 0 : weight[0] = 0.0e0;
389 : 0 : double maxdist = -1e100;
390 : 0 : double totweight = 0.0e0;
391 [ # # ]: 0 : for(j=0; j<num_close; j++) {
392 [ # # ][ # # ]: 0 : cp = close_points[j]->coordinates();
393 [ # # ]: 0 : point->transform_to_local( cp, coords[j] );
394 [ # # ][ # # ]: 0 : weight[j] = sqrt( FacetEvalToolUtils::sqr(coords[j].x()) + FacetEvalToolUtils::sqr(coords[j].y()) );
[ # # ][ # # ]
395 [ # # ]: 0 : if (weight[j] > maxdist) maxdist = weight[j];
396 : : }
397 : 0 : maxdist *= 1.1e0;
398 [ # # ]: 0 : for (j=0; j<num_close; j++) {
399 [ # # ]: 0 : weight[j] = FacetEvalToolUtils::sqr((maxdist-weight[j])/(maxdist*weight[j]));
400 : 0 : totweight += weight[j];
401 : : }
402 : :
403 : : // fill up the matrices
404 : :
405 : : //lhs.set_to_identity();
406 : : //rhs.set_to_identity();
407 : : //coef.set_to_identity();
408 : :
409 : : double dx, dy, wjdx, wjdy, dx2, dy2, dxdy, dz;
410 [ # # ]: 0 : for (j=0; j<num_close; j++) {
411 : 0 : weight[j] /= totweight;
412 : 0 : weight[j] = 1; // BUGFIX: ignore weights for now and reset weights to 1
413 [ # # ]: 0 : dx = /*-*/ coords[j].x(); //BUGFIX: Why we need -ve coords?
414 [ # # ]: 0 : dy = /*-*/ coords[j].y(); //BUGFIX: Why we need -ve coords?
415 : 0 : wjdx = weight[j] * dx;
416 : 0 : wjdy = weight[j] * dy;
417 [ # # ]: 0 : dx2 = FacetEvalToolUtils::sqr( dx );
418 [ # # ]: 0 : dy2 = FacetEvalToolUtils::sqr( dy );
419 : 0 : dxdy = dx * dy;
420 [ # # ]: 0 : dz = coords[j].z();
421 : :
422 [ # # ]: 0 : lhs.add( 0, 0, wjdx * dx );
423 [ # # ]: 0 : lhs.add( 0, 1, wjdx * dy );
424 [ # # ]: 0 : lhs.add( 0, 2, wjdx * dx2 );
425 [ # # ]: 0 : lhs.add( 0, 3, wjdx * dxdy );
426 [ # # ]: 0 : lhs.add( 0, 4, wjdx * dy2 );
427 [ # # ]: 0 : rhs.add( 0, 0, wjdx * dz ); // BUGFIX: dz was missing
428 : :
429 [ # # ]: 0 : lhs.add( 1, 1, wjdy * dy );
430 [ # # ]: 0 : lhs.add( 1, 2, wjdy * dx2 );
431 [ # # ]: 0 : lhs.add( 1, 3, wjdy * dxdy );
432 [ # # ]: 0 : lhs.add( 1, 4, wjdy * dx * dy2 );
433 [ # # ]: 0 : rhs.add( 1, 0, wjdy * dz ); // BUGFIX: dz was missing
434 : :
435 [ # # ]: 0 : lhs.add( 2, 2, wjdx * dx2 * dx );
436 [ # # ]: 0 : lhs.add( 2, 3, wjdx * dx2 * dy );
437 [ # # ]: 0 : lhs.add( 2, 4, wjdx * dx * dy2 );
438 [ # # ]: 0 : rhs.add( 2, 0, wjdx * dx * dz );// BUGFIX: dz was missing
439 : :
440 [ # # ]: 0 : lhs.add( 3, 3, wjdx * dx * dy2 );
441 [ # # ]: 0 : lhs.add( 3, 4, wjdx * dy * dy2 );
442 [ # # ]: 0 : rhs.add( 3, 0, wjdx * dy * dz ); // BUGFIX: dz was missing
443 : :
444 [ # # ]: 0 : lhs.add( 4, 4, wjdy * dy * dy2 );
445 [ # # ]: 0 : rhs.add( 4, 0, wjdy * dy * dz ); // BUGFIX: dz was missing
446 : : }
447 [ # # ][ # # ]: 0 : lhs.set( 1, 0, lhs.get(0,1) );
448 [ # # ][ # # ]: 0 : lhs.set( 2, 0, lhs.get(0,2) );
449 [ # # ][ # # ]: 0 : lhs.set( 2, 1, lhs.get(1,2) );
450 [ # # ][ # # ]: 0 : lhs.set( 3, 0, lhs.get(0,3) );
451 [ # # ][ # # ]: 0 : lhs.set( 3, 1, lhs.get(1,3) );
452 [ # # ][ # # ]: 0 : lhs.set( 3, 2, lhs.get(2,3) );
453 [ # # ][ # # ]: 0 : lhs.set( 4, 0, lhs.get(0,4) );
454 [ # # ][ # # ]: 0 : lhs.set( 4, 1, lhs.get(1,4) );
455 [ # # ][ # # ]: 0 : lhs.set( 4, 2, lhs.get(2,4) );
456 [ # # ][ # # ]: 0 : lhs.set( 4, 3, lhs.get(3,4) );
457 : :
458 : : // solve the system
459 : :
460 [ # # ]: 0 : status = lhs.solveNxN( rhs, coef );
461 [ # # ]: 0 : if (status != CUBIT_SUCCESS) {
462 : 0 : return status;
463 : : }
464 : :
465 : : // store the quadric coefficents with the point
466 : :
467 [ # # ]: 0 : point->coef_vector( coef );
468 : : }
469 [ # # ]: 0 : return CUBIT_SUCCESS;
470 : : }
471 : :
472 : : //===========================================================================
473 : : //Function Name: get_close_points
474 : : //
475 : : //Member Type: PRIVATE
476 : : //Descriptoin: get a list of points close to the current point on the facets
477 : : // return a minimum of min_close and a maximum of max_close points
478 : : //===========================================================================
479 : 0 : CubitStatus FacetEvalTool::get_close_points(
480 : : CubitPoint *point,
481 : : CubitPoint **close_point,
482 : : int &num_close,
483 : : int max_close,
484 : : int min_close )
485 : : {
486 : :
487 : : // get the points immediately adjacent
488 : :
489 [ # # ]: 0 : DLIList<CubitFacet*> adj_facet_list;
490 [ # # ]: 0 : point->facets(adj_facet_list);
491 [ # # ][ # # ]: 0 : if (adj_facet_list.size() > max_close) {
492 : 0 : return CUBIT_FAILURE;
493 : : }
494 [ # # ]: 0 : point->adjacent_points(close_point, num_close);
495 : :
496 : : // if we don't have enough yet, then go to the next level
497 : :
498 [ # # ]: 0 : if (num_close < min_close) {
499 : : CubitPoint *cpoint[100];
500 : : CubitPoint *cur_point;
501 : : CubitBoolean found;
502 : : int nclose, cnclose;
503 : 0 : nclose = num_close;
504 [ # # ]: 0 : for(int i=0; i<num_close; i++) {
505 : 0 : cur_point = close_point[i];
506 [ # # ]: 0 : DLIList<CubitFacet*> cadj_facet_list;
507 [ # # ]: 0 : cur_point->facets(cadj_facet_list);
508 [ # # ][ # # ]: 0 : if (cadj_facet_list.size() + nclose > max_close) {
509 : 0 : return CUBIT_FAILURE;
510 : : }
511 [ # # ]: 0 : cur_point->adjacent_points(cpoint,cnclose);
512 [ # # ][ # # ]: 0 : for(int k=0; k<cnclose; k++) {
[ # # ]
513 : :
514 : : // check that it is not already on the list
515 : :
516 : 0 : found = CUBIT_FALSE;
517 [ # # ][ # # ]: 0 : for(int l=0; l<nclose && !found; l++) {
518 [ # # ][ # # ]: 0 : if (close_point[l] == cpoint[k] ||
519 : 0 : cpoint[k] == point) {
520 : 0 : found = CUBIT_TRUE;
521 : : }
522 : : }
523 [ # # ]: 0 : if (!found) {
524 : :
525 : : // make sure the normal at this point is not more than 90 degrees
526 : : // from the normal at the point
527 : :
528 [ # # ][ # # ]: 0 : double dot = cpoint[k]->normal() % point->normal();
[ # # ]
529 [ # # ]: 0 : if (dot > GEOMETRY_RESABS) {
530 : :
531 : : // add the point to the list
532 : :
533 : 0 : close_point[nclose++] = cpoint[k];
534 : : }
535 : : }
536 : : }
537 : 0 : }
538 : 0 : num_close = nclose;
539 [ # # ]: 0 : if (num_close < min_close) {
540 : 0 : return CUBIT_FAILURE;
541 : : }
542 : : }
543 [ # # ]: 0 : return CUBIT_SUCCESS;
544 : : }
545 : :
546 : : //===========================================================================
547 : : //Function Name: mark edge pairs
548 : : //
549 : : //Member Type: PRIVATE
550 : : //Descriptoin: For the given point, we loop over the attached boundary
551 : : // edges and figure out which pairs should be C1 continuous.
552 : : // A given edge is paired with at most two other edges. The
553 : : // other edges (or the other edge) is stored on a tool data
554 : : // to be used later.
555 : : //===========================================================================
556 : 264 : CubitStatus FacetEvalTool::mark_edge_pairs(CubitPoint* point)
557 : : {
558 : : int i,j;
559 [ + - ]: 264 : DLIList<CubitFacetEdge*> edge_list;
560 [ + - ][ + - ]: 528 : DLIList<CubitFacetEdge*> edge_list_init;
561 : 264 : CubitPoint* prev_point = NULL;
562 : 264 : CubitPoint* next_point = NULL;
563 : 264 : CubitFacetEdge* prev_edge = NULL;
564 : 264 : CubitFacetEdge* next_edge = NULL;
565 [ + - ][ + - ]: 264 : CubitVector e0, e1;
566 : : double current_dot;
567 [ + - ]: 264 : point->edges(edge_list_init);
568 : 264 : TDFacetBoundaryEdge *td_fbe = NULL;
569 : :
570 : : // make a list with only the boundary edges
571 [ + - ][ + + ]: 1452 : for(i=0; i< edge_list_init.size(); i++){
572 [ + - ]: 1188 : prev_edge = edge_list_init.get_and_step();
573 [ + - ]: 1188 : td_fbe = TDFacetBoundaryEdge::get_facet_boundary_edge(prev_edge);
574 [ + + ]: 1188 : if(td_fbe)
575 [ + - ]: 792 : edge_list.append(prev_edge);
576 : : }
577 : 264 : prev_edge=NULL;
578 : : //if there is one or none, we don't need to bother looking
579 : : // for pairs.
580 [ + - ][ - + ]: 264 : if(edge_list.size() < 2)
581 : 0 : return CUBIT_SUCCESS;
582 [ + - ]: 264 : int num_edges=edge_list.size();
583 [ + - ][ + - ]: 528 : std::vector<int> other_index(num_edges);
584 [ + - ][ + - ]: 528 : std::vector<double> dot_array(num_edges);
585 : 264 : bool pair_exists = false;
586 : : // loop over the edges. For each edge, find the edge that would
587 : : // make the best other edge to pair this one with
588 [ + + ]: 1056 : for(i = 0; i<num_edges; i++){
589 [ + - ]: 792 : other_index[i]=-1;
590 [ + - ]: 792 : dot_array[i]=-1.0;
591 : :
592 [ + - ]: 792 : prev_edge=edge_list[i];
593 [ + - ]: 792 : prev_point = prev_edge->other_point(point);
594 : : // now figure out which other edge would be the best pair with
595 : : // this one. This could be sped up by only looking at the edges
596 : : // i+1 to num_edges, but ever this more exhaustive search is not
597 : : // guaranteed to make an optimal pairing... we take the longer
598 : : // approach hoping it will give slightly better results for
599 : : // hard problems...
600 [ + + ]: 3168 : for(j = 0; j<num_edges; j++){
601 [ + + ]: 2376 : if(j!=i){
602 [ + - ]: 1584 : next_edge = edge_list[j];
603 [ + - ]: 1584 : next_point = next_edge->other_point(point);
604 [ + - ][ + - ]: 1584 : e0 = point->coordinates() - next_point->coordinates();
[ + - ][ + - ]
605 [ + - ][ + - ]: 1584 : e1 = prev_point->coordinates() - point->coordinates();
[ + - ][ + - ]
606 [ + - ]: 1584 : e0.normalize();
607 [ + - ]: 1584 : e1.normalize();
608 [ + - ]: 1584 : current_dot = e0%e1;
609 : : //if the current dot satisfies the feature angle criterion
610 : : // and is better than any other we've seen so far for this
611 : : // given edge, save it.
612 [ - + ][ # # ]: 1584 : if(current_dot >= minDot && current_dot > dot_array[i]){
[ # # ][ - + ]
613 [ # # ]: 0 : dot_array[i]=current_dot;
614 [ # # ]: 0 : other_index[i]=j;
615 : 0 : pair_exists=true;//keep track of whether a pair has been saved
616 : : }
617 : : }
618 : : }
619 : : }
620 : : //if there aren't any pairs, don't bother moving forward.
621 [ + - ]: 264 : if(!pair_exists)
622 : 264 : return CUBIT_SUCCESS;
623 : : //now find the best pair. That is the pair with biggest
624 : : // dot product. Then find the next, and so on until we
625 : : // are done.
626 : 0 : double best_this_time = CUBIT_DBL_MAX;
627 : 0 : int best_index=-1;
628 : : // given num_edges > 2 and each edge is paired with at most
629 : : // one other edge at this node, there can't be more than
630 : : // num_edges - 1 pairs. Actually, num_edges / 2, but
631 : : // it is just a safety check anyway, so num_edges-1 will work.
632 [ # # ][ # # ]: 0 : for(i=0;i<num_edges-1 && best_this_time >= minDot; i++){
633 : 0 : best_this_time = -1.0;
634 : 0 : best_index = -1;
635 : : //loop over and find the biggest dot
636 [ # # ]: 0 : for(j=0;j<num_edges; j++){
637 [ # # ][ # # ]: 0 : if(dot_array[j] > best_this_time){
638 [ # # ]: 0 : best_this_time = dot_array[j];
639 : 0 : best_index = j;
640 : : }
641 : : }
642 : : //if we found a pair that we can make C1
643 [ # # ]: 0 : if(best_index >= 0){
644 : : //Don't let the above loop find either of these again
645 : : // (unless they are the 'other' in the pair).
646 [ # # ]: 0 : dot_array[best_index] = -1.0;
647 [ # # ][ # # ]: 0 : dot_array[other_index[best_index]] = -1.0;
648 : :
649 : : //First, make sure the other in the pair hasn't already
650 : : // been used.
651 [ # # ]: 0 : CubitFacetEdge* edge_1 = edge_list[best_index];
652 [ # # ][ # # ]: 0 : CubitFacetEdge* edge_2 = edge_list[other_index[best_index]];
653 [ # # ]: 0 : td_fbe = TDFacetBoundaryEdge::get_facet_boundary_edge( edge_2 );
654 [ # # ]: 0 : if(!td_fbe){
655 [ # # ][ # # ]: 0 : PRINT_ERROR("Expected a tool data.\n");
[ # # ][ # # ]
656 : 0 : return CUBIT_FAILURE;
657 : : }
658 : : //figure out whether it should be stored as prev or next
659 [ # # ][ # # ]: 0 : if(point == edge_2->point(0)){
660 : : //if something has already been stored here, then
661 : : // need to skip this pair
662 [ # # ][ # # ]: 0 : if(td_fbe->get_prev_edge())
663 : 0 : edge_1 = NULL;
664 : : else
665 [ # # ]: 0 : td_fbe->set_prev_edge(edge_1);
666 : : }
667 : : else{
668 [ # # ][ # # ]: 0 : if(td_fbe->get_next_edge())
669 : 0 : edge_1 = NULL;
670 : : else
671 [ # # ]: 0 : td_fbe->set_next_edge(edge_1);
672 : : }
673 : : //edge_1 will be NULL if we decided above to skip this pair
674 [ # # ]: 0 : if(edge_1){//otherwise save edge_2 in the appropriate spot
675 : : // for edge_1's tool data.
676 [ # # ]: 0 : td_fbe = TDFacetBoundaryEdge::get_facet_boundary_edge( edge_1 );
677 [ # # ]: 0 : if(!td_fbe){
678 [ # # ][ # # ]: 0 : PRINT_ERROR("Expected another tool data.\n");
[ # # ][ # # ]
679 : 0 : return CUBIT_FAILURE;
680 : : }
681 : :
682 [ # # ][ # # ]: 0 : if(point == edge_1->point(0))
683 [ # # ]: 0 : td_fbe->set_prev_edge(edge_2);
684 : : else
685 [ # # ]: 0 : td_fbe->set_next_edge(edge_2);
686 : : }
687 : : }
688 : :
689 : : }
690 [ + - ]: 264 : return CUBIT_SUCCESS;
691 : : }
692 : :
693 : :
694 : :
695 : : //===========================================================================
696 : : //Function Name: pair_edges
697 : : //
698 : : //Member Type: PRIVATE
699 : : //Descriptoin: Basically, loops over the points and calls
700 : : // mark_edge_pairs
701 : : //===========================================================================
702 : 66 : CubitStatus FacetEvalTool::pair_edges()
703 : : {
704 : 66 : CubitStatus status = CUBIT_SUCCESS;
705 : : int i;
706 : 66 : CubitPoint* point = NULL;
707 [ + + ]: 330 : for( i=0;i<myPointList.size(); i++){
708 : 264 : point = myPointList.get_and_step();
709 : 264 : status = mark_edge_pairs(point);
710 [ - + ]: 264 : if(status!=CUBIT_SUCCESS)
711 : 0 : return status;
712 : : }
713 : 66 : return status;
714 : : }
715 : :
716 : : //===========================================================================
717 : : //Function Name: init_bezier_surface
718 : : //
719 : : //Member Type: PRIVATE
720 : : //Descriptoin: compute the surface control points
721 : : //===========================================================================
722 : 66 : CubitStatus FacetEvalTool::init_bezier_surface()
723 : : {
724 : : // initialize the edges
725 : :
726 : : int i;
727 : 66 : CubitStatus status = CUBIT_SUCCESS;
728 : : CubitFacetEdge *edge;
729 : : // figure out which edges should be paired for C1 continuity.
730 : 66 : status = pair_edges();
731 [ - + ]: 66 : if(status!=CUBIT_SUCCESS)
732 : 0 : return status;
733 : :
734 [ + + ][ + - ]: 396 : for (i=0; i<myEdgeList.size() && status == CUBIT_SUCCESS; i++) {
[ + + ]
735 : 330 : edge = myEdgeList.get_and_step();
736 : 330 : status = init_bezier_edge( edge, minDot );
737 [ - + ]: 330 : if (status != CUBIT_SUCCESS)
738 : 0 : return status;
739 : : }
740 : 66 : int mydebug = 0;
741 [ - + ]: 66 : if (mydebug)
742 : : {
743 : 0 : debug_draw_bezier_edges();
744 : 0 : dview();
745 : : }
746 : :
747 : : // initialize the facets
748 : :
749 [ + - ]: 66 : if (status == CUBIT_SUCCESS) {
750 : : CubitFacet *facet;
751 [ + + ][ + - ]: 198 : for (i=0; i<myFacetList.size() && status == CUBIT_SUCCESS; i++) {
[ + + ]
752 : 132 : facet = myFacetList.get_and_step();
753 : 132 : status = init_bezier_facet( facet );
754 : : }
755 : : }
756 [ - + ]: 66 : if(status != CUBIT_SUCCESS){
757 [ # # ][ # # ]: 0 : PRINT_ERROR("Problem initializing bezier facet.\n");
758 : 0 : return status;
759 : : }
760 : :
761 : : // reset the bounding box to account for new control points
762 : :
763 : 66 : reset_bounding_box();
764 : :
765 : : //draw_eval_bezier_facet( myFacetList.get_and_step() );
766 : : // draw_bezier_facet( myFacetList.get_and_step() );
767 : :
768 : 66 : return status;
769 : : }
770 : :
771 : : //===========================================================================
772 : : //Function Name: init_bezier_edge
773 : : //
774 : : //Member Type: PRIVATE
775 : : //Descriptoin: compute the control points for an edge
776 : : //===========================================================================
777 : 330 : CubitStatus FacetEvalTool::init_bezier_edge( CubitFacetEdge *edge,
778 : : double min_dot )
779 : : {
780 : 330 : CubitStatus stat = CUBIT_SUCCESS;
781 : :
782 : : TDFacetBoundaryEdge *td_bfe =
783 : 330 : TDFacetBoundaryEdge::get_facet_boundary_edge( edge );
784 [ + + ]: 330 : if (td_bfe)
785 : : {
786 : 264 : stat = td_bfe->init_control_points( min_dot );
787 [ - + ]: 264 : if (stat != CUBIT_SUCCESS)
788 : 0 : return stat;
789 : 264 : stat = td_bfe->merge_control_points();
790 : : }
791 : : else
792 : : {
793 [ + - ][ + + ]: 264 : CubitVector ctrl_pts[3];
794 [ + - ][ + - ]: 66 : CubitVector P0 = edge->point(0)->coordinates();
795 [ + - ][ + - ]: 66 : CubitVector P3 = edge->point(1)->coordinates();
796 [ + - ][ + - ]: 66 : CubitVector N0 = edge->point(0)->normal( edge );
797 [ + - ][ + - ]: 66 : CubitVector N3 = edge->point(1)->normal( edge );
798 [ + - ][ + - ]: 66 : CubitVector T0, T3;
799 [ + - ][ - + ]: 66 : if (edge->num_adj_facets() <= 1)
800 : : {
801 [ # # ]: 0 : stat = compute_curve_tangent( edge, min_dot, T0, T3 );
802 [ # # ]: 0 : if (stat != CUBIT_SUCCESS)
803 : 0 : return stat;
804 : : }
805 : : else
806 : : {
807 [ + - ][ + - ]: 66 : T0 = P3 - P0;
808 [ + - ]: 66 : T0.normalize();
809 [ + - ]: 66 : T3 = T0;
810 : : }
811 [ + - ]: 66 : stat = init_edge_control_points( P0, P3, N0, N3, T0, T3, ctrl_pts );
812 [ - + ]: 66 : if (stat != CUBIT_SUCCESS)
813 : 0 : return stat;
814 [ + - ]: 66 : edge->control_points( ctrl_pts, 4 );
815 : : }
816 : 330 : return stat;
817 : : }
818 : :
819 : :
820 : : //===========================================================================
821 : : //Function Name: init_edge_control_points
822 : : //
823 : : //Member Type: PRIVATE
824 : : //Descriptoin: compute the control points for an edge
825 : : //===========================================================================
826 : 594 : CubitStatus FacetEvalTool::init_edge_control_points( CubitVector &P0,
827 : : CubitVector &P3,
828 : : CubitVector &N0,
829 : : CubitVector &N3,
830 : : CubitVector &T0,
831 : : CubitVector &T3,
832 : : CubitVector Pi[3] )
833 : : {
834 [ + - ][ + + ]: 2970 : CubitVector Vi[4];
835 [ + - ]: 594 : Vi[0] = P0;
836 [ + - ]: 594 : Vi[3] = P3;
837 [ + - ]: 594 : double di = P0.distance_between( P3 );
838 [ + - ]: 594 : double ai = N0 % N3;
839 [ + - ]: 594 : double ai0 = N0 % T0;
840 [ + - ]: 594 : double ai3 = N3 % T3;
841 [ + - ]: 594 : double denom = 4 - FacetEvalToolUtils::sqr(ai);
842 [ - + ]: 594 : if (fabs(denom) < 1e-20) {
843 : 0 : return CUBIT_FAILURE;
844 : : }
845 : 594 : double row = 6.0e0 * (2.0e0 * ai0 + ai * ai3) / denom;
846 : 594 : double omega = 6.0e0 * (2.0e0 * ai3 + ai * ai0) / denom;
847 [ + - ][ + - ]: 594 : Vi[1] = Vi[0] + (di * (((6.0e0 * T0) - ((2.0e0 * row) * N0) + (omega * N3)) / 18.0e0));
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
848 [ + - ][ + - ]: 594 : Vi[2] = Vi[3] - (di * (((6.0e0 * T3) + (row * N0) - ((2.0e0 * omega) * N3)) / 18.0e0));
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
849 [ + - ][ + + ]: 2376 : CubitVector Wi[3];
850 [ + - ][ + - ]: 594 : Wi[0] = Vi[1] - Vi[0];
851 [ + - ][ + - ]: 594 : Wi[1] = Vi[2] - Vi[1];
852 [ + - ][ + - ]: 594 : Wi[2] = Vi[3] - Vi[2];
853 : :
854 [ + - ][ + - ]: 594 : Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
[ + - ][ + - ]
855 [ + - ][ + - ]: 594 : Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
[ + - ][ + - ]
856 [ + - ][ + - ]: 594 : Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
[ + - ][ + - ]
857 : :
858 : 594 : return CUBIT_SUCCESS;
859 : : }
860 : :
861 : : //===========================================================================
862 : : //Function Name: init_edge_control_points_single
863 : : //
864 : : //Member Type: PRIVATE
865 : : //Descriptoin: compute the control points for an edge without normals
866 : : //===========================================================================
867 : 0 : CubitStatus FacetEvalTool::init_edge_control_points_single(
868 : : CubitVector &P0,
869 : : CubitVector &P3,
870 : : CubitVector &T0,
871 : : CubitVector &T3,
872 : : CubitVector Pi[3] )
873 : : {
874 [ # # ][ # # ]: 0 : CubitVector Vi[4];
875 [ # # ]: 0 : Vi[0] = P0;
876 [ # # ]: 0 : Vi[3] = P3;
877 [ # # ]: 0 : double di = P0.distance_between( P3 );
878 [ # # ]: 0 : double ai = T0 % T3;
879 [ # # ]: 0 : double denom = 4 - FacetEvalToolUtils::sqr(ai);
880 [ # # ]: 0 : if (fabs(denom) < 1e-20) {
881 : 0 : return CUBIT_FAILURE;
882 : : }
883 [ # # ][ # # ]: 0 : Vi[1] = Vi[0] + (di * (((6.0e0 * T0)) / 18.0e0));
[ # # ][ # # ]
[ # # ]
884 [ # # ][ # # ]: 0 : Vi[2] = Vi[3] - (di * (((6.0e0 * T3)) / 18.0e0));
[ # # ][ # # ]
[ # # ]
885 : :
886 [ # # ][ # # ]: 0 : Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
[ # # ][ # # ]
887 [ # # ][ # # ]: 0 : Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
[ # # ][ # # ]
888 [ # # ][ # # ]: 0 : Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
[ # # ][ # # ]
889 : :
890 : 0 : return CUBIT_SUCCESS;
891 : : }
892 : :
893 : : //===========================================================================
894 : : //Function Name: init_bezier_facet
895 : : //
896 : : //Member Type: PRIVATE
897 : : //Descriptoin: compute the control points for a facet
898 : : //===========================================================================
899 : 264 : CubitStatus FacetEvalTool::init_bezier_facet( CubitFacet *facet )
900 : : {
901 : :
902 : 264 : CubitStatus stat = CUBIT_SUCCESS;
903 [ + - ][ + + ]: 5016 : CubitVector P[3][5];
[ + + ]
904 [ + - ][ + + ]: 3432 : CubitVector N[6], G[6];
[ + - ][ + + ]
905 [ + - ]: 264 : stat = facet->get_edge_control_points( P );
906 [ - + ]: 264 : if (stat != CUBIT_SUCCESS)
907 : 0 : return stat;
908 : :
909 : : // retreive the normals from the edge points. Note we duplicate the
910 : : // pointer data here only because of the edge_use
911 : :
912 : 264 : int mydebug = 0;
913 [ - + ]: 264 : if (mydebug)
914 : : {
915 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
916 [ # # ]: 0 : dfdraw(facet);
917 [ # # ]: 0 : dview();
918 : : }
919 [ + + ]: 1056 : for (int i=0; i<3; i++) {
920 [ + - ]: 792 : CubitFacetEdge *edge = facet->edge(i);
921 [ + - ][ + + ]: 792 : if (facet->edge_use(i) == 1) {
922 [ + - ][ + - ]: 396 : N[i*2] = edge->point(0)->normal( facet );
[ + - ]
923 [ + - ][ + - ]: 396 : N[i*2+1] = edge->point(1)->normal( facet );
[ + - ]
924 : : }
925 : : else {
926 [ + - ][ + - ]: 396 : N[i*2] = edge->point(1)->normal( facet );
[ + - ]
927 [ + - ][ + - ]: 396 : N[i*2+1] = edge->point(0)->normal( facet );
[ + - ]
928 : : }
929 : : }
930 : :
931 : : // init the facet control points.
932 : :
933 [ + - ]: 264 : stat = init_facet_control_points( N, P, G );
934 [ - + ]: 264 : if (stat != CUBIT_SUCCESS)
935 : 0 : return stat;
936 [ + - ]: 264 : facet->set_control_points ( G );
937 [ + - ]: 264 : facet->update_bezier_bounding_box();
938 : 264 : return stat;
939 : : }
940 : :
941 : : //===========================================================================
942 : : //Function Name: init_facet_control_points
943 : : //
944 : : //Member Type: PRIVATE
945 : : //Descriptoin: compute the control points for a facet
946 : : //===========================================================================
947 : 264 : CubitStatus FacetEvalTool::init_facet_control_points(
948 : : CubitVector N[6], // vertex normals (per edge)
949 : : CubitVector P[3][5], // edge control points
950 : : CubitVector G[6] ) // return internal control points
951 : : {
952 [ + - ][ + + ]: 3960 : CubitVector Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ][ + + ]
953 : : double denom;
954 : : double lambda[2], mu[2];
955 : :
956 : 264 : CubitStatus stat = CUBIT_SUCCESS;
957 : :
958 [ + + ]: 1056 : for (int i=0; i<3; i++) {
959 [ + - ]: 792 : N0 = N[i*2];
960 [ + - ]: 792 : N3 = N[i*2+1];
961 [ + - ]: 792 : Vi[0] = P[i][0];
962 [ + - ][ + - ]: 792 : Vi[1] = (P[i][1] - 0.25 * P[i][0]) / 0.75;
[ + - ][ + - ]
963 [ + - ][ + - ]: 792 : Vi[2] = (P[i][3] - 0.25 * P[i][4]) / 0.75;
[ + - ][ + - ]
964 [ + - ]: 792 : Vi[3] = P[i][4];
965 [ + - ][ + - ]: 792 : Wi[0] = Vi[1] - Vi[0];
966 [ + - ][ + - ]: 792 : Wi[1] = Vi[2] - Vi[1];
967 [ + - ][ + - ]: 792 : Wi[2] = Vi[3] - Vi[2];
968 [ + - ][ + - ]: 792 : Di[0] = P[(i+2)%3][3] - 0.5*(P[i][1] + P[i][0]);
[ + - ][ + - ]
969 [ + - ][ + - ]: 792 : Di[3] = P[(i+1)%3][1] - 0.5*(P[i][4] + P[i][3]);
[ + - ][ + - ]
970 [ + - ][ - + ]: 792 : if(Wi[0].length() == 0.0)
971 : : {
972 [ # # ][ # # ]: 0 : Ai[0] = N0 * Wi[0];
973 [ # # ]: 0 : lambda[0] = Di[0] % Wi[0];
974 : : }
975 : : else
976 : : {
977 [ + - ][ + - ]: 792 : Ai[0] = (N0 * Wi[0]) / Wi[0].length();
[ + - ][ + - ]
978 [ + - ][ + - ]: 792 : lambda[0] = (Di[0] % Wi[0]) / (Wi[0] % Wi[0]);
979 : : }
980 : :
981 [ + - ][ - + ]: 792 : if(Wi[2].length() == 0.0)
982 : : {
983 [ # # ][ # # ]: 0 : Ai[2] = N3 * Wi[2];
984 [ # # ]: 0 : lambda[1] = Di[3] % Wi[2];
985 : : }
986 : : else
987 : : {
988 [ + - ][ + - ]: 792 : Ai[2] = (N3 * Wi[2]) / Wi[2].length();
[ + - ][ + - ]
989 [ + - ][ + - ]: 792 : lambda[1] = (Di[3] % Wi[2]) / (Wi[2] % Wi[2]);
990 : : }
991 [ + - ][ + - ]: 792 : Ai[1] = Ai[0] + Ai[2];
992 [ + - ]: 792 : denom = Ai[1].length();
993 [ + - ]: 792 : if (denom > 0)
994 [ + - ]: 792 : Ai[1] /= denom;
995 [ + - ]: 792 : mu[0] = (Di[0] % Ai[0]);
996 [ + - ]: 792 : mu[1] = (Di[3] % Ai[2]);
997 [ + - ][ + - ]: 792 : G[i*2] = 0.5 * (P[i][1] + P[i][2]) +
[ + - ]
998 [ + - ][ + - ]: 1584 : 0.66666666666666 * lambda[0] * Wi[1] +
999 [ + - ][ + - ]: 1584 : 0.33333333333333 * lambda[1] * Wi[0] +
1000 [ + - ][ + - ]: 1584 : 0.66666666666666 * mu[0] * Ai[1] +
1001 [ + - ][ + - ]: 1584 : 0.33333333333333 * mu[1] * Ai[0];
1002 [ + - ][ + - ]: 792 : G[i*2+1] = 0.5 * (P[i][2] + P[i][3]) +
[ + - ]
1003 [ + - ][ + - ]: 1584 : 0.33333333333333 * lambda[0] * Wi[2] +
1004 [ + - ][ + - ]: 1584 : 0.66666666666666 * lambda[1] * Wi[1] +
1005 [ + - ][ + - ]: 1584 : 0.33333333333333 * mu[0] * Ai[2] +
1006 [ + - ][ + - ]: 1584 : 0.66666666666666 * mu[1] * Ai[1];
1007 : : }
1008 : 264 : return stat;
1009 : : }
1010 : :
1011 : : //===========================================================================
1012 : : //Function Name: eval_bezier_patch
1013 : : //
1014 : : //Member Type: PRIVATE
1015 : : //Descriptoin: evaluate the bezier patch defined at a facet
1016 : : //===========================================================================
1017 : 0 : CubitStatus FacetEvalTool::eval_bezier_patch( CubitFacet *facet,
1018 : : CubitVector &areacoord,
1019 : : CubitVector &pt)
1020 : : {
1021 : : // interpolate internal control points
1022 : :
1023 [ # # ][ # # ]: 0 : CubitVector gctrl_pts[6];
1024 [ # # ]: 0 : facet->get_control_points( gctrl_pts );
1025 [ # # ][ # # ]: 0 : CubitVector P_facet[3];
1026 [ # # ][ # # ]: 0 : if (fabs(areacoord.y() + areacoord.z()) < 1.0e-6) {
[ # # ]
1027 [ # # ][ # # ]: 0 : pt = facet->point(0)->coordinates();
[ # # ]
1028 : 0 : return CUBIT_SUCCESS;
1029 : : }
1030 [ # # ][ # # ]: 0 : if (fabs(areacoord.x() + areacoord.z()) < 1.0e-6) {
[ # # ]
1031 [ # # ][ # # ]: 0 : pt = facet->point(1)->coordinates();
[ # # ]
1032 : 0 : return CUBIT_SUCCESS;
1033 : : }
1034 [ # # ][ # # ]: 0 : if (fabs(areacoord.x() + areacoord.y()) < 1.0e-6) {
[ # # ]
1035 [ # # ][ # # ]: 0 : pt = facet->point(2)->coordinates();
[ # # ]
1036 : 0 : return CUBIT_SUCCESS;
1037 : : }
1038 : :
1039 : : //2,1,1
1040 [ # # ][ # # ]: 0 : P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
[ # # ]
1041 [ # # ][ # # ]: 0 : (areacoord.y() * gctrl_pts[3] +
[ # # ]
1042 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[4]);
[ # # ]
1043 : : //1,2,1
1044 [ # # ][ # # ]: 0 : P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
[ # # ]
1045 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[0] +
[ # # ]
1046 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[5]);
[ # # ]
1047 : : //1,1,2
1048 [ # # ][ # # ]: 0 : P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
[ # # ]
1049 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[1] +
[ # # ]
1050 [ # # ][ # # ]: 0 : areacoord.y() * gctrl_pts[2]);
[ # # ]
1051 : :
1052 : : // sum the contribution from each of the control points
1053 : :
1054 [ # # ]: 0 : pt.set(0.0e0, 0.0e0, 0.0e0);
1055 : : CubitFacetEdge *edge;
1056 [ # # ]: 0 : edge = facet->edge( 2 );
1057 [ # # ][ # # ]: 0 : CubitVector ctrl_pts[5];
1058 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1059 : :
1060 : : //i=4; j=0; k=0;
1061 [ # # ][ # # ]: 0 : double B = FacetEvalToolUtils::quart(areacoord.x());
1062 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[0];
1063 : :
1064 : : //i=3; j=1; k=0;
1065 [ # # ][ # # ]: 0 : B = 4.0 * FacetEvalToolUtils::cube(areacoord.x()) * areacoord.y();
[ # # ]
1066 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[1];
1067 : :
1068 : : //i=2; j=2; k=0;
1069 [ # # ][ # # ]: 0 : B = 6.0 * FacetEvalToolUtils::sqr(areacoord.x()) * FacetEvalToolUtils::sqr(areacoord.y());
[ # # ][ # # ]
1070 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[2];
1071 : :
1072 : : //i=1; j=3; k=0;
1073 [ # # ][ # # ]: 0 : B = 4.0 * areacoord.x() * FacetEvalToolUtils::cube(areacoord.y());
[ # # ]
1074 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[3];
1075 : :
1076 [ # # ]: 0 : edge = facet->edge( 0 );
1077 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1078 : :
1079 : : //i=0; j=4; k=0;
1080 [ # # ][ # # ]: 0 : B = FacetEvalToolUtils::quart(areacoord.y());
1081 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[0];
1082 : :
1083 : : //i=0; j=3; k=1;
1084 [ # # ][ # # ]: 0 : B = 4.0 * FacetEvalToolUtils::cube(areacoord.y()) * areacoord.z();
[ # # ]
1085 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[1];
1086 : :
1087 : : //i=0; j=2; k=2;
1088 [ # # ][ # # ]: 0 : B = 6.0 * FacetEvalToolUtils::sqr(areacoord.y()) * FacetEvalToolUtils::sqr(areacoord.z());
[ # # ][ # # ]
1089 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[2];
1090 : :
1091 : : //i=0; j=1; k=3;
1092 [ # # ][ # # ]: 0 : B = 4.0 * areacoord.y() * FacetEvalToolUtils::cube(areacoord.z());
[ # # ]
1093 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[3];
1094 : :
1095 [ # # ]: 0 : edge = facet->edge( 1 );
1096 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1097 : :
1098 : : //i=0; j=0; k=4;
1099 [ # # ][ # # ]: 0 : B = FacetEvalToolUtils::quart(areacoord.z());
1100 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[0];
1101 : :
1102 : : //i=1; j=0; k=3;
1103 [ # # ][ # # ]: 0 : B = 4.0 * areacoord.x() * FacetEvalToolUtils::cube(areacoord.z());
[ # # ]
1104 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[1];
1105 : :
1106 : : //i=2; j=0; k=2;
1107 [ # # ][ # # ]: 0 : B = 6.0 * FacetEvalToolUtils::sqr(areacoord.x()) * FacetEvalToolUtils::sqr(areacoord.z());
[ # # ][ # # ]
1108 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[2];
1109 : :
1110 : : //i=3; j=0; k=1;
1111 [ # # ][ # # ]: 0 : B = 4.0 * FacetEvalToolUtils::cube(areacoord.x()) * areacoord.z();
[ # # ]
1112 [ # # ][ # # ]: 0 : pt += B * ctrl_pts[3];
1113 : :
1114 : : //i=2; j=1; k=1;
1115 [ # # ][ # # ]: 0 : B = 12.0 * FacetEvalToolUtils::sqr(areacoord.x()) * areacoord.y() * areacoord.z();
[ # # ][ # # ]
1116 [ # # ][ # # ]: 0 : pt += B * P_facet[0];
1117 : :
1118 : : //i=1; j=2; k=1;
1119 [ # # ][ # # ]: 0 : B = 12.0 * areacoord.x() * FacetEvalToolUtils::sqr(areacoord.y()) * areacoord.z();
[ # # ][ # # ]
1120 [ # # ][ # # ]: 0 : pt += B * P_facet[1];
1121 : :
1122 : : //i=1; j=1; k=2;
1123 [ # # ][ # # ]: 0 : B = 12.0 * areacoord.x() * areacoord.y() * FacetEvalToolUtils::sqr(areacoord.z());
[ # # ][ # # ]
1124 [ # # ][ # # ]: 0 : pt += B * P_facet[2];
1125 : :
1126 : 0 : return CUBIT_SUCCESS;
1127 : : }
1128 : :
1129 : : //===========================================================================
1130 : : //Function Name: eval_bezier_patch_normal
1131 : : //
1132 : : //Member Type: PRIVATE
1133 : : //Descriptoin: evaluate the bezier patch defined at a facet
1134 : : //===========================================================================
1135 : 0 : CubitStatus FacetEvalTool::eval_bezier_patch_normal( CubitFacet *facet,
1136 : : CubitVector &areacoord,
1137 : : CubitVector &normal )
1138 : : {
1139 : : // interpolate internal control points
1140 : :
1141 [ # # ][ # # ]: 0 : CubitVector gctrl_pts[6];
1142 [ # # ]: 0 : facet->get_control_points( gctrl_pts );
1143 [ # # ][ # # ]: 0 : if (fabs(areacoord.y() + areacoord.z()) < 1.0e-6) {
[ # # ]
1144 [ # # ][ # # ]: 0 : normal = facet->point(0)->normal(facet);
[ # # ]
1145 : 0 : return CUBIT_SUCCESS;
1146 : : }
1147 [ # # ][ # # ]: 0 : if (fabs(areacoord.x() + areacoord.z()) < 1.0e-6) {
[ # # ]
1148 [ # # ][ # # ]: 0 : normal = facet->point(1)->normal(facet);
[ # # ]
1149 : 0 : return CUBIT_SUCCESS;
1150 : : }
1151 [ # # ][ # # ]: 0 : if (fabs(areacoord.x() + areacoord.y()) < 1.0e-6) {
[ # # ]
1152 [ # # ][ # # ]: 0 : normal = facet->point(2)->normal(facet);
[ # # ]
1153 : 0 : return CUBIT_SUCCESS;
1154 : : }
1155 : :
1156 : : // compute the hodograph of the quartic Gregory patch
1157 : :
1158 [ # # ][ # # ]: 0 : CubitVector Nijk[10];
1159 [ # # ]: 0 : hodograph(facet,areacoord,Nijk);
1160 : :
1161 : : // sum the contribution from each of the control points
1162 : :
1163 [ # # ]: 0 : normal.set(0.0e0, 0.0e0, 0.0e0);
1164 : :
1165 : : //i=3; j=0; k=0;
1166 : 0 : double Bsum = 0.0;
1167 [ # # ][ # # ]: 0 : double B = FacetEvalToolUtils::cube(areacoord.x());
1168 : 0 : Bsum += B;
1169 [ # # ][ # # ]: 0 : normal += B * Nijk[0];
1170 : :
1171 : : //i=2; j=1; k=0;
1172 [ # # ][ # # ]: 0 : B = 3.0 * FacetEvalToolUtils::sqr(areacoord.x()) * areacoord.y();
[ # # ]
1173 : 0 : Bsum += B;
1174 [ # # ][ # # ]: 0 : normal += B * Nijk[1];
1175 : :
1176 : : //i=1; j=2; k=0;
1177 [ # # ][ # # ]: 0 : B = 3.0 * areacoord.x() * FacetEvalToolUtils::sqr(areacoord.y());
[ # # ]
1178 : 0 : Bsum += B;
1179 [ # # ][ # # ]: 0 : normal += B * Nijk[2];
1180 : :
1181 : : //i=0; j=3; k=0;
1182 [ # # ][ # # ]: 0 : B = FacetEvalToolUtils::cube(areacoord.y());
1183 : 0 : Bsum += B;
1184 [ # # ][ # # ]: 0 : normal += B * Nijk[3];
1185 : :
1186 : : //i=2; j=0; k=1;
1187 [ # # ][ # # ]: 0 : B = 3.0 * FacetEvalToolUtils::sqr(areacoord.x()) * areacoord.z();
[ # # ]
1188 : 0 : Bsum += B;
1189 [ # # ][ # # ]: 0 : normal += B * Nijk[4];
1190 : :
1191 : : //i=1; j=1; k=1;
1192 [ # # ][ # # ]: 0 : B = 6.0 * areacoord.x() * areacoord.y() * areacoord.z();
[ # # ]
1193 : 0 : Bsum += B;
1194 [ # # ][ # # ]: 0 : normal += B * Nijk[5];
1195 : :
1196 : : //i=0; j=2; k=1;
1197 [ # # ][ # # ]: 0 : B = 3.0 * FacetEvalToolUtils::sqr(areacoord.y()) * areacoord.z();
[ # # ]
1198 : 0 : Bsum += B;
1199 [ # # ][ # # ]: 0 : normal += B * Nijk[6];
1200 : :
1201 : : //i=1; j=0; k=2;
1202 [ # # ][ # # ]: 0 : B = 3.0 * areacoord.x() * FacetEvalToolUtils::sqr(areacoord.z());
[ # # ]
1203 : 0 : Bsum += B;
1204 [ # # ][ # # ]: 0 : normal += B * Nijk[7];
1205 : :
1206 : : //i=0; j=1; k=2;
1207 [ # # ][ # # ]: 0 : B = 3.0 * areacoord.y() * FacetEvalToolUtils::sqr(areacoord.z());
[ # # ]
1208 : 0 : Bsum += B;
1209 [ # # ][ # # ]: 0 : normal += B * Nijk[8];
1210 : :
1211 : : //i=0; j=0; k=3;
1212 [ # # ][ # # ]: 0 : B = FacetEvalToolUtils::cube(areacoord.z());
1213 : 0 : Bsum += B;
1214 [ # # ][ # # ]: 0 : normal += B * Nijk[9];
1215 : :
1216 [ # # ]: 0 : assert(fabs(Bsum - 1.0) < 1e-9);
1217 : :
1218 [ # # ]: 0 : normal.normalize();
1219 : :
1220 : 0 : return CUBIT_SUCCESS;
1221 : : }
1222 : :
1223 : : //===========================================================================
1224 : : //Function Name: hodograph
1225 : : //
1226 : : //Member Type: PUBLIC
1227 : : //Description: get the hodograph control points for the facet
1228 : : //Note: This is a triangle cubic patch that is defined by the
1229 : : // normals of quartic facet control point lattice. Returned coordinates
1230 : : // in Nijk are defined by the following diagram
1231 : : //
1232 : : //
1233 : : // *9 index polar
1234 : : // / \ 0 300 point(0)
1235 : : // / \ 1 210
1236 : : // 7*-----*8 2 120
1237 : : // / \ / \ 3 030 point(1)
1238 : : // / \ / \ 4 201
1239 : : // 4*----5*-----*6 5 111
1240 : : // / \ / \ / \ 6 021
1241 : : // / \ / \ / \ 7 102
1242 : : // *-----*-----*-----* 8 012
1243 : : // 0 1 2 3 9 003 point(2)
1244 : : //
1245 : : //===========================================================================
1246 : 0 : CubitStatus FacetEvalTool::hodograph( CubitFacet *facet,
1247 : : CubitVector &areacoord,
1248 : : CubitVector Nijk[10] )
1249 : : {
1250 : :
1251 : : // compute the control points on the interior of the patch based on areacoord
1252 : :
1253 [ # # ][ # # ]: 0 : CubitVector gctrl_pts[6];
1254 [ # # ]: 0 : facet->get_control_points( gctrl_pts );
1255 [ # # ][ # # ]: 0 : CubitVector P_facet[3];
1256 : :
1257 : : //2,1,1
1258 [ # # ][ # # ]: 0 : P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
[ # # ]
1259 [ # # ][ # # ]: 0 : (areacoord.y() * gctrl_pts[3] +
[ # # ]
1260 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[4]);
[ # # ]
1261 : : //1,2,1
1262 [ # # ][ # # ]: 0 : P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
[ # # ]
1263 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[0] +
[ # # ]
1264 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[5]);
[ # # ]
1265 : : //1,1,2
1266 [ # # ][ # # ]: 0 : P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
[ # # ]
1267 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[1] +
[ # # ]
1268 [ # # ][ # # ]: 0 : areacoord.y() * gctrl_pts[2]);
[ # # ]
1269 : :
1270 : : // corner control points are just the normals at the points
1271 : :
1272 : : //3, 0, 0
1273 [ # # ][ # # ]: 0 : Nijk[0] = facet->point(0)->normal(facet);
[ # # ]
1274 : : //0, 3, 0
1275 [ # # ][ # # ]: 0 : Nijk[3] = facet->point(1)->normal(facet);
[ # # ]
1276 : : //0, 0, 3
1277 [ # # ][ # # ]: 0 : Nijk[9] = facet->point(2)->normal(facet);
[ # # ]
1278 : :
1279 : : // fill in the boundary control points. Define as the normal to the local
1280 : : // triangle formed by the quartic control point lattice
1281 : :
1282 : : CubitFacetEdge *edge;
1283 [ # # ]: 0 : edge = facet->edge( 2 );
1284 [ # # ][ # # ]: 0 : CubitVector ctrl_pts[5];
1285 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1286 : :
1287 : : //2, 1, 0
1288 [ # # ][ # # ]: 0 : Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
[ # # ][ # # ]
1289 [ # # ]: 0 : Nijk[1].normalize();
1290 : :
1291 : : //1, 2, 0
1292 [ # # ][ # # ]: 0 : Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
[ # # ][ # # ]
1293 [ # # ]: 0 : Nijk[2].normalize();
1294 : :
1295 [ # # ]: 0 : edge = facet->edge( 0 );
1296 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1297 : :
1298 : : //0, 2, 1
1299 [ # # ][ # # ]: 0 : Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
[ # # ][ # # ]
1300 [ # # ]: 0 : Nijk[6].normalize();
1301 : :
1302 : : //0, 1, 2
1303 [ # # ][ # # ]: 0 : Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
[ # # ][ # # ]
1304 [ # # ]: 0 : Nijk[8].normalize();
1305 : :
1306 [ # # ]: 0 : edge = facet->edge( 1 );
1307 [ # # ]: 0 : edge->control_points(facet, ctrl_pts);
1308 : :
1309 : : //1, 0, 2
1310 [ # # ][ # # ]: 0 : Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
[ # # ][ # # ]
1311 [ # # ]: 0 : Nijk[7].normalize();
1312 : :
1313 : : //2, 0, 1
1314 [ # # ][ # # ]: 0 : Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
[ # # ][ # # ]
1315 [ # # ]: 0 : Nijk[4].normalize();
1316 : :
1317 : : //1, 1, 1
1318 [ # # ][ # # ]: 0 : Nijk[5] = (P_facet[1] - P_facet[0]) * (P_facet[2] - P_facet[0]);
[ # # ][ # # ]
1319 [ # # ]: 0 : Nijk[5].normalize();
1320 : :
1321 : 0 : int mydebug = 0;
1322 [ # # ]: 0 : if (mydebug)
1323 : : {
1324 : : #define ONE_THIRD 0.33333333333333333333333333333333333
1325 : : #define TWO_THIRDS .666666666666666666666666666666666666667
1326 : : int ii;
1327 [ # # ][ # # ]: 0 : CubitVector ac, pt;
1328 [ # # ]: 0 : for(ii=0; ii<10; ii++)
1329 : : {
1330 [ # # # # : 0 : switch(ii)
# # # # #
# # ]
1331 : : {
1332 [ # # ]: 0 : case 0: ac.set(1.0, 0.0, 0.0 ); break;
1333 [ # # ]: 0 : case 1: ac.set(TWO_THIRDS, ONE_THIRD, 0.0 ); break;
1334 [ # # ]: 0 : case 2: ac.set(ONE_THIRD, TWO_THIRDS, 0.0 ); break;
1335 [ # # ]: 0 : case 3: ac.set(0.0, 1.0, 0.0 ); break;
1336 [ # # ]: 0 : case 4: ac.set(TWO_THIRDS, 0.0, ONE_THIRD ); break;
1337 [ # # ]: 0 : case 5: ac.set(ONE_THIRD, ONE_THIRD, ONE_THIRD ); break;
1338 [ # # ]: 0 : case 6: ac.set(0.0, TWO_THIRDS, ONE_THIRD ); break;
1339 [ # # ]: 0 : case 7: ac.set(ONE_THIRD, 0.0, TWO_THIRDS); break;
1340 [ # # ]: 0 : case 8: ac.set(0.0, ONE_THIRD, TWO_THIRDS); break;
1341 [ # # ]: 0 : case 9: ac.set(0.0, 0.0, 1.0 ); break;
1342 : : }
1343 [ # # ]: 0 : eval_bezier_patch(facet, ac, pt);
1344 [ # # ]: 0 : dray(pt,Nijk[ii],1.0);
1345 : : }
1346 [ # # ]: 0 : dview();
1347 : : }
1348 : :
1349 : 0 : return CUBIT_SUCCESS;
1350 : : }
1351 : :
1352 : :
1353 : : //===========================================================================
1354 : : //Function Name: project_to_patch
1355 : : //
1356 : : //Member Type: PUBLIC
1357 : : //Descriptoin: Project a point to a bezier patch. Pass in the area
1358 : : // of the point projected to the linear facet. Function
1359 : : // assumes that the point is contained within the patch -
1360 : : // if not, it will project to one of its edges.
1361 : : //===========================================================================
1362 : 0 : CubitStatus FacetEvalTool::project_to_patch(
1363 : : CubitFacet *facet, // (IN) the facet where the patch is defined
1364 : : CubitVector &ac, // (IN) area coordinate initial guess (from linear facet)
1365 : : const CubitVector &pt, // (IN) point we are projecting to patch
1366 : : CubitVector &eval_pt, // (OUT) The projected point
1367 : : CubitVector *eval_norm, // (OUT) normal at evaluated point
1368 : : CubitBoolean &outside, // (OUT) the closest point on patch to pt is on an edge
1369 : : double compare_tol, // (IN) comparison tolerance
1370 : : int edge_id ) // (IN) only used if this is to be projected to one
1371 : : // of the edges. Otherwise, should be -1
1372 : : {
1373 : 0 : CubitStatus status = CUBIT_SUCCESS;
1374 : :
1375 : : // see if we are at a vertex
1376 : :
1377 : : #define INCR 0.01
1378 : 0 : const double tol = compare_tol;
1379 : :
1380 [ # # ][ # # ]: 0 : if (is_at_vertex( facet, pt, ac, compare_tol, eval_pt, eval_norm ))
1381 : : {
1382 : 0 : outside = CUBIT_FALSE;
1383 : 0 : return CUBIT_SUCCESS;
1384 : : }
1385 : :
1386 : : // check if the start ac is inside the patch -if not, then move it there
1387 : :
1388 : 0 : int nout = 0;
1389 : 0 : const double atol = 0.001;
1390 [ # # ][ # # ]: 0 : if(move_ac_inside( ac, atol ))
1391 : 0 : nout++;
1392 : :
1393 : 0 : int diverge = 0;
1394 : 0 : int iter = 0;
1395 [ # # ]: 0 : CubitVector newpt;
1396 [ # # ]: 0 : eval_bezier_patch(facet, ac, newpt);
1397 [ # # ]: 0 : CubitVector move = pt - newpt;
1398 [ # # ]: 0 : double lastdist = move.length();
1399 : 0 : double bestdist = lastdist;
1400 [ # # ]: 0 : CubitVector bestac = ac;
1401 [ # # ]: 0 : CubitVector bestpt = newpt;
1402 [ # # ]: 0 : CubitVector bestnorm;
1403 : :
1404 : : // If we are already close enough, then return now
1405 : :
1406 [ # # ][ # # ]: 0 : if (lastdist <= tol && !eval_norm && nout == 0) {
[ # # ]
1407 [ # # ]: 0 : eval_pt = pt;
1408 : 0 : outside = CUBIT_FALSE;
1409 : 0 : return status;
1410 : : }
1411 : :
1412 : : double ratio, mag, umove, vmove, det, distnew, movedist;
1413 [ # # ]: 0 : CubitVector lastpt = newpt;
1414 [ # # ]: 0 : CubitVector lastac = ac;
1415 [ # # ]: 0 : CubitVector norm;
1416 [ # # ][ # # ]: 0 : CubitVector xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1417 [ # # ][ # # ]: 0 : CubitVector du, dv, newac;
[ # # ]
1418 : 0 : CubitBoolean done = CUBIT_FALSE;
1419 [ # # ]: 0 : while (!done) {
1420 : :
1421 : : // We will be locating the projected point within the u,v,w coordinate
1422 : : // system of the triangular bezier patch. Since u+v+w=1, the components
1423 : : // are not linearly independant. We will choose only two of the
1424 : : // coordinates to use and compute the third.
1425 : :
1426 : : int system;
1427 [ # # ][ # # ]: 0 : if (lastac.x() >= lastac.y() && lastac.x() >= lastac.z()) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1428 : 0 : system = 0;
1429 : : }
1430 [ # # ][ # # ]: 0 : else if (lastac.y() >= lastac.z()) {
[ # # ]
1431 : 0 : system = 1;
1432 : : }
1433 : : else {
1434 : 0 : system = 2;
1435 : : }
1436 : :
1437 : : // compute the surface derivatives with respect to each
1438 : : // of the barycentric coordinates
1439 : :
1440 : :
1441 [ # # ][ # # ]: 0 : if (system == 1 || system == 2) {
1442 [ # # ][ # # ]: 0 : xac.x( lastac.x() + INCR );
1443 [ # # ][ # # ]: 0 : if (lastac.y() + lastac.z() == 0.0)
[ # # ]
1444 : 0 : return CUBIT_FAILURE;
1445 [ # # ][ # # ]: 0 : ratio = lastac.z() / (lastac.y() + lastac.z());
[ # # ]
1446 [ # # ][ # # ]: 0 : xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
1447 [ # # ][ # # ]: 0 : xac.z( 1.0 - xac.x() - xac.y() );
[ # # ]
1448 [ # # ]: 0 : eval_bezier_patch(facet, xac, xpt);
1449 [ # # ][ # # ]: 0 : xvec = xpt - lastpt;
1450 [ # # ]: 0 : xvec /= INCR;
1451 : : }
1452 [ # # ][ # # ]: 0 : if (system == 0 || system == 2) {
1453 [ # # ][ # # ]: 0 : yac.y( lastac.y() + INCR );
1454 [ # # ][ # # ]: 0 : if (lastac.x() + lastac.z() == 0.0)
[ # # ]
1455 : 0 : return CUBIT_FAILURE;
1456 [ # # ][ # # ]: 0 : ratio = lastac.z() / (lastac.x() + lastac.z());
[ # # ]
1457 [ # # ][ # # ]: 0 : yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
1458 [ # # ][ # # ]: 0 : yac.z( 1.0 - yac.x() - yac.y() );
[ # # ]
1459 [ # # ]: 0 : eval_bezier_patch(facet, yac, ypt);
1460 [ # # ][ # # ]: 0 : yvec = ypt - lastpt;
1461 [ # # ]: 0 : yvec /= INCR;
1462 : : }
1463 [ # # ][ # # ]: 0 : if (system == 0 || system == 1) {
1464 [ # # ][ # # ]: 0 : zac.z( lastac.z() + INCR );
1465 [ # # ][ # # ]: 0 : if (lastac.x() + lastac.y() == 0.0)
[ # # ]
1466 : 0 : return CUBIT_FAILURE;
1467 [ # # ][ # # ]: 0 : ratio = lastac.y() / (lastac.x() + lastac.y());
[ # # ]
1468 [ # # ][ # # ]: 0 : zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
1469 [ # # ][ # # ]: 0 : zac.y( 1.0 - zac.x() - zac.z() );
[ # # ]
1470 [ # # ]: 0 : eval_bezier_patch(facet, zac, zpt);
1471 [ # # ][ # # ]: 0 : zvec = zpt - lastpt;
1472 [ # # ]: 0 : zvec /= INCR;
1473 : : }
1474 : :
1475 : : // compute the surface normal
1476 : :
1477 [ # # # # ]: 0 : switch (system) {
1478 : : case 0:
1479 [ # # ]: 0 : du = yvec;
1480 [ # # ]: 0 : dv = zvec;
1481 : 0 : break;
1482 : : case 1:
1483 [ # # ]: 0 : du = zvec;
1484 [ # # ]: 0 : dv = xvec;
1485 : 0 : break;
1486 : : case 2:
1487 [ # # ]: 0 : du = xvec;
1488 [ # # ]: 0 : dv = yvec;
1489 : 0 : break;
1490 : : }
1491 [ # # ][ # # ]: 0 : norm = du * dv;
1492 [ # # ]: 0 : mag = norm.length();
1493 [ # # ]: 0 : if (mag < DBL_EPSILON) {
1494 : 0 : return CUBIT_FAILURE;
1495 : : // do something else here (it is likely a flat triangle -
1496 : : // so try evaluating just an edge of the bezier patch)
1497 : : }
1498 [ # # ]: 0 : norm /= mag;
1499 [ # # ]: 0 : if (iter == 0)
1500 [ # # ]: 0 : bestnorm = norm;
1501 : :
1502 : : // project the move vector to the tangent plane
1503 : :
1504 [ # # ][ # # ]: 0 : move = (norm * move) * norm;
[ # # ]
1505 : :
1506 : : // compute an equivalent u-v-w vector
1507 : :
1508 [ # # ][ # # ]: 0 : CubitVector absnorm( fabs(norm.x()), fabs(norm.y()), fabs(norm.z()) );
[ # # ][ # # ]
1509 [ # # ][ # # ]: 0 : if (absnorm.z() >= absnorm.y() && absnorm.z() >= absnorm.x()) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1510 [ # # ][ # # ]: 0 : det = du.x() * dv.y() - dv.x() * du.y();
[ # # ][ # # ]
1511 [ # # ]: 0 : if (fabs(det) <= DBL_EPSILON) {
1512 : 0 : return CUBIT_FAILURE; // do something else here
1513 : : }
1514 [ # # ][ # # ]: 0 : umove = (move.x() * dv.y() - dv.x() * move.y()) / det;
[ # # ][ # # ]
1515 [ # # ][ # # ]: 0 : vmove = (du.x() * move.y() - move.x() * du.y()) / det;
[ # # ][ # # ]
1516 : : }
1517 [ # # ][ # # ]: 0 : else if (absnorm.y() >= absnorm.z() && absnorm.y() >= absnorm.x()) {
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1518 [ # # ][ # # ]: 0 : det = du.x() * dv.z() - dv.x() * du.z();
[ # # ][ # # ]
1519 [ # # ]: 0 : if (fabs(det) <= DBL_EPSILON) {
1520 : 0 : return CUBIT_FAILURE;
1521 : : }
1522 [ # # ][ # # ]: 0 : umove = (move.x() * dv.z() - dv.x() * move.z()) / det;
[ # # ][ # # ]
1523 [ # # ][ # # ]: 0 : vmove = (du.x() * move.z() - move.x() * du.z()) / det;
[ # # ][ # # ]
1524 : : }
1525 : : else {
1526 [ # # ][ # # ]: 0 : det = du.y() * dv.z() - dv.y() * du.z();
[ # # ][ # # ]
1527 [ # # ]: 0 : if (fabs(det) <= DBL_EPSILON) {
1528 : 0 : return CUBIT_FAILURE;
1529 : : }
1530 [ # # ][ # # ]: 0 : umove = (move.y() * dv.z() - dv.y() * move.z()) / det;
[ # # ][ # # ]
1531 [ # # ][ # # ]: 0 : vmove = (du.y() * move.z() - move.y() * du.z()) / det;
[ # # ][ # # ]
1532 : : }
1533 : :
1534 : : /* === compute the new u-v coords and evaluate surface at new location */
1535 : :
1536 [ # # # # ]: 0 : switch (system) {
1537 : : case 0:
1538 [ # # ][ # # ]: 0 : newac.y( lastac.y() + umove );
1539 [ # # ][ # # ]: 0 : newac.z( lastac.z() + vmove );
1540 [ # # ][ # # ]: 0 : newac.x( 1.0 - newac.y() - newac.z() );
[ # # ]
1541 : 0 : break;
1542 : : case 1:
1543 [ # # ][ # # ]: 0 : newac.z( lastac.z() + umove );
1544 [ # # ][ # # ]: 0 : newac.x( lastac.x() + vmove );
1545 [ # # ][ # # ]: 0 : newac.y( 1.0 - newac.z() - newac.x() );
[ # # ]
1546 : 0 : break;
1547 : : case 2:
1548 [ # # ][ # # ]: 0 : newac.x( lastac.x() + umove );
1549 [ # # ][ # # ]: 0 : newac.y( lastac.y() + vmove );
1550 [ # # ][ # # ]: 0 : newac.z( 1.0 - newac.x() - newac.y() );
[ # # ]
1551 : 0 : break;
1552 : : }
1553 : :
1554 : : // Keep it inside the patch
1555 : :
1556 [ # # ][ # # ]: 0 : if ( newac.x() >= -atol &&
[ # # ]
1557 [ # # ][ # # ]: 0 : newac.y() >= -atol &&
[ # # ]
1558 [ # # ]: 0 : newac.z() >= -atol) {
1559 : 0 : nout = 0;
1560 : : }
1561 : : else {
1562 [ # # ][ # # ]: 0 : if (move_ac_inside( newac, atol ) == CUBIT_TRUE)
1563 : 0 : nout++;
1564 : : }
1565 : :
1566 : : // Evaluate at the new location
1567 : :
1568 [ # # ]: 0 : if (edge_id != -1)
1569 [ # # ]: 0 : ac_at_edge( newac, newac, edge_id ); // move to edge first
1570 [ # # ]: 0 : eval_bezier_patch(facet, newac, newpt);
1571 : :
1572 : : // Check for convergence
1573 : :
1574 [ # # ]: 0 : distnew = pt.distance_between(newpt);
1575 [ # # ][ # # ]: 0 : move = newpt - lastpt;
1576 [ # # ]: 0 : movedist = move.length();
1577 [ # # ][ # # ]: 0 : if (movedist < tol ||
1578 : : distnew < tol ) {
1579 : 0 : done = CUBIT_TRUE;
1580 [ # # ]: 0 : if (distnew < bestdist)
1581 : : {
1582 : 0 : bestdist = distnew;
1583 [ # # ]: 0 : bestac = newac;
1584 [ # # ]: 0 : bestpt = newpt;
1585 [ # # ]: 0 : bestnorm = norm;
1586 : : }
1587 : : }
1588 : : else {
1589 : :
1590 : : // don't allow more than 30 iterations
1591 : :
1592 : 0 : iter++;
1593 [ # # ]: 0 : if (iter > 30) {
1594 : : //if (movedist > tol * 100.0) nout=1;
1595 : 0 : done = CUBIT_TRUE;
1596 : : }
1597 : :
1598 : : // Check for divergence - don't allow more than 5 divergent
1599 : : // iterations
1600 : :
1601 [ # # ]: 0 : if (distnew > lastdist) {
1602 : 0 : diverge++;
1603 [ # # ]: 0 : if (diverge > 10) {
1604 : 0 : done = CUBIT_TRUE;
1605 : : //if (movedist > tol * 100.0) nout=1;
1606 : : }
1607 : : }
1608 : :
1609 : : // Check if we are continuing to project outside the facet.
1610 : : // If so, then stop now
1611 : :
1612 [ # # ]: 0 : if (nout > 3) {
1613 : 0 : done = CUBIT_TRUE;
1614 : : }
1615 : :
1616 : : // set up for next iteration
1617 : :
1618 [ # # ]: 0 : if (!done) {
1619 [ # # ]: 0 : if (distnew < bestdist)
1620 : : {
1621 : 0 : bestdist = distnew;
1622 [ # # ]: 0 : bestac = newac;
1623 [ # # ]: 0 : bestpt = newpt;
1624 [ # # ]: 0 : bestnorm = norm;
1625 : : }
1626 : 0 : lastdist = distnew;
1627 [ # # ]: 0 : lastpt = newpt;
1628 [ # # ]: 0 : lastac = newac;
1629 [ # # ][ # # ]: 0 : move = pt - lastpt;
1630 : : }
1631 : : }
1632 : : }
1633 : :
1634 [ # # ]: 0 : eval_pt = bestpt;
1635 [ # # ]: 0 : if (eval_norm) {
1636 [ # # ]: 0 : *eval_norm = bestnorm;
1637 : : }
1638 [ # # ]: 0 : outside = (nout > 0) ? CUBIT_TRUE : CUBIT_FALSE;
1639 [ # # ]: 0 : ac = bestac;
1640 : :
1641 : 0 : return status;
1642 : : }
1643 : :
1644 : : //===========================================================================
1645 : : //Function Name: is_at_vertex
1646 : : //
1647 : : //Member Type: PRIVATE
1648 : : //Description: determine if the point is at one of the facet's vertices
1649 : : //===========================================================================
1650 : 0 : CubitBoolean FacetEvalTool::is_at_vertex(
1651 : : CubitFacet *facet, // (IN) facet we are evaluating
1652 : : const CubitVector &pt, // (IN) the point
1653 : : CubitVector &ac, // (IN) the ac of the point on the facet plane
1654 : : double compare_tol, // (IN) return TRUE of closer than this
1655 : : CubitVector &eval_pt, // (OUT) location at vertex if TRUE
1656 : : CubitVector *eval_norm_ptr ) // (OUT) normal at vertex if TRUE
1657 : : {
1658 : : double dist;
1659 [ # # ]: 0 : CubitVector vert_loc;
1660 : 0 : const double actol = 0.1;
1661 [ # # ][ # # ]: 0 : if (fabs(ac.x()) < actol && fabs(ac.y()) < actol) {
[ # # ][ # # ]
[ # # ]
1662 [ # # ][ # # ]: 0 : vert_loc = facet->point(2)->coordinates();
[ # # ]
1663 [ # # ]: 0 : dist = pt.distance_between( vert_loc );
1664 [ # # ]: 0 : if (dist <= compare_tol)
1665 : : {
1666 [ # # ]: 0 : eval_pt = vert_loc;
1667 [ # # ]: 0 : if (eval_norm_ptr)
1668 : : {
1669 [ # # ][ # # ]: 0 : *eval_norm_ptr = facet->point(2)->normal( facet );
[ # # ]
1670 : : }
1671 : 0 : return CUBIT_TRUE;
1672 : : }
1673 : : }
1674 : :
1675 [ # # ][ # # ]: 0 : if (fabs(ac.x()) < actol && fabs(ac.z()) < actol) {
[ # # ][ # # ]
[ # # ]
1676 [ # # ][ # # ]: 0 : vert_loc = facet->point(1)->coordinates();
[ # # ]
1677 [ # # ]: 0 : dist = pt.distance_between( vert_loc );
1678 [ # # ]: 0 : if (dist <= compare_tol)
1679 : : {
1680 [ # # ]: 0 : eval_pt = vert_loc;
1681 [ # # ]: 0 : if (eval_norm_ptr)
1682 : : {
1683 [ # # ][ # # ]: 0 : *eval_norm_ptr = facet->point(1)->normal( facet );
[ # # ]
1684 : : }
1685 : 0 : return CUBIT_TRUE;
1686 : : }
1687 : : }
1688 : :
1689 [ # # ][ # # ]: 0 : if (fabs(ac.y()) < actol && fabs(ac.z()) < actol) {
[ # # ][ # # ]
[ # # ]
1690 [ # # ][ # # ]: 0 : vert_loc = facet->point(0)->coordinates();
[ # # ]
1691 [ # # ]: 0 : dist = pt.distance_between( vert_loc );
1692 [ # # ]: 0 : if (dist <= compare_tol)
1693 : : {
1694 [ # # ]: 0 : eval_pt = vert_loc;
1695 [ # # ]: 0 : if (eval_norm_ptr)
1696 : : {
1697 [ # # ][ # # ]: 0 : *eval_norm_ptr = facet->point(0)->normal( facet );
[ # # ]
1698 : : }
1699 : 0 : return CUBIT_TRUE;
1700 : : }
1701 : : }
1702 : :
1703 : 0 : return CUBIT_FALSE;
1704 : : }
1705 : :
1706 : : //===========================================================================
1707 : : //Function Name: ac_at_edge
1708 : : //
1709 : : //Member Type: PRIVATE
1710 : : //Description: determine the area coordinate of the facet at the edge
1711 : : //===========================================================================
1712 : 0 : void FacetEvalTool::ac_at_edge( CubitVector &fac, // facet area coordinate
1713 : : CubitVector &eac, // edge area coordinate
1714 : : int edge_id ) // id of edge
1715 : : {
1716 : 0 : double u = 0.0, v = 0.0, w = 0.0;
1717 [ # # # # ]: 0 : switch (edge_id)
1718 : : {
1719 : : case 0:
1720 : 0 : u = 0.0;
1721 : 0 : v = fac.y() / (fac.y() + fac.z());
1722 : 0 : w = 1.0 - v;
1723 : 0 : break;
1724 : : case 1:
1725 : 0 : u = fac.x() / (fac.x() + fac.z());
1726 : 0 : v = 0.0;
1727 : 0 : w = 1.0 - u;
1728 : 0 : break;
1729 : : case 2:
1730 : 0 : u = fac.x() / (fac.x() + fac.y());
1731 : 0 : v = 1.0 - u;
1732 : 0 : w = 0.0;
1733 : 0 : break;
1734 : : default:
1735 : 0 : assert(0);
1736 : : break;
1737 : : }
1738 : 0 : eac.set(u, v, w);
1739 : 0 : }
1740 : :
1741 : : //===========================================================================
1742 : : //Function Name: move_ac_inside
1743 : : //
1744 : : //Member Type: PRIVATE
1745 : : //Description: find the closest area coordinate to the boundary of the
1746 : : // patch if any of its components are < 0
1747 : : // Return if the ac was modified.
1748 : : //===========================================================================
1749 : 0 : CubitBoolean FacetEvalTool::move_ac_inside( CubitVector &ac, double tol )
1750 : : {
1751 : 0 : int nout = 0;
1752 [ # # ]: 0 : if (ac.x() < -tol) {
1753 : 0 : ac.x( 0.0 );
1754 : 0 : ac.y( ac.y() / (ac.y() + ac.z()) );
1755 : 0 : ac.z( 1.0 - ac.y() );
1756 : 0 : nout++;
1757 : : }
1758 [ # # ]: 0 : if (ac.y() < -tol) {
1759 : 0 : ac.y( 0.0 );
1760 : 0 : ac.x( ac.x() / (ac.x() + ac.z()) );
1761 : 0 : ac.z( 1.0 - ac.x() );
1762 : 0 : nout++;
1763 : : }
1764 [ # # ]: 0 : if (ac.z() < -tol) {
1765 : 0 : ac.z( 0.0 );
1766 : 0 : ac.x( ac.x() / (ac.x() + ac.y()) );
1767 : 0 : ac.y( 1.0 - ac.x() );
1768 : 0 : nout++;
1769 : : }
1770 [ # # ]: 0 : return (nout > 0) ? CUBIT_TRUE : CUBIT_FALSE;
1771 : : }
1772 : :
1773 : 473 : bool FacetEvalTool::have_data_to_calculate_bbox(void)
1774 : : {
1775 [ - + ][ # # ]: 473 : if(myPointList.size() > 0 ||
[ + - ]
1776 [ # # ]: 0 : (interpOrder == 4 &&
1777 [ # # ]: 0 : (myEdgeList.size() > 0 ||
1778 : 0 : myFacetList.size() > 0)))
1779 : : {
1780 : 473 : return true;
1781 : : }
1782 : 0 : return false;
1783 : : }
1784 : :
1785 : : //===========================================================================
1786 : : //Function Name: reset_bounding_box
1787 : : //
1788 : : //Member Type: PUBLIC
1789 : : //Description: Calculates the bounding box of the surface (also sets
1790 : : // the compareTol and grid search).
1791 : : //===========================================================================
1792 : 473 : void FacetEvalTool::reset_bounding_box()
1793 : : {
1794 [ + - ]: 473 : if(have_data_to_calculate_bbox())
1795 : : {
1796 [ + + ]: 473 : if (myBBox != NULL)
1797 : : {
1798 [ + - ]: 132 : delete myBBox;
1799 : 132 : myBBox = NULL;
1800 : : }
1801 : :
1802 : 473 : bounding_box();
1803 : :
1804 : 473 : double diag = sqrt(FacetEvalToolUtils::sqr(myBBox->x_range()) +
1805 : 473 : FacetEvalToolUtils::sqr(myBBox->y_range()) +
1806 : 473 : FacetEvalToolUtils::sqr(myBBox->z_range()));
1807 : 473 : compareTol = 1.0e-3 * diag;
1808 : :
1809 : 473 : set_up_grid_search( diag );
1810 : : }
1811 : 473 : }
1812 : :
1813 : :
1814 : : //===========================================================================
1815 : : //Function Name: bounding_box
1816 : : //
1817 : : //Member Type: PUBLIC
1818 : : //Descriptoin: Calculates the bounding box of the surface.
1819 : : //===========================================================================
1820 : 2838 : CubitBox FacetEvalTool::bounding_box()
1821 : : {
1822 [ + + ]: 2838 : if ( !myBBox )
1823 : : {
1824 : : int ii;
1825 : : CubitPoint *point_ptr;
1826 : : double x, y, z;
1827 : 605 : double x_min = CUBIT_DBL_MAX, x_max = -CUBIT_DBL_MAX;
1828 : 605 : double y_min = CUBIT_DBL_MAX, y_max = -CUBIT_DBL_MAX;
1829 : 605 : double z_min = CUBIT_DBL_MAX, z_max = -CUBIT_DBL_MAX;
1830 [ + - ][ + + ]: 3377 : for ( ii = myPointList.size(); ii > 0; ii-- )
1831 : : {
1832 [ + - ]: 2772 : point_ptr = myPointList.get_and_step();
1833 [ + - ]: 2772 : x = point_ptr->x();
1834 [ + - ]: 2772 : y = point_ptr->y();
1835 [ + - ]: 2772 : z = point_ptr->z();
1836 [ + + ]: 2772 : if ( x < x_min )
1837 : 825 : x_min = x;
1838 [ + + ]: 2772 : if ( y < y_min )
1839 : 715 : y_min = y;
1840 [ + + ]: 2772 : if ( z < z_min )
1841 : 759 : z_min = z;
1842 [ + + ]: 2772 : if ( x > x_max )
1843 : 880 : x_max = x;
1844 [ + + ]: 2772 : if ( y > y_max )
1845 : 1056 : y_max = y;
1846 [ + + ]: 2772 : if ( z > z_max )
1847 : 1067 : z_max = z;
1848 : : }
1849 [ + + ]: 605 : if (interpOrder == 4)
1850 : : {
1851 : : CubitFacetEdge *edge_ptr;
1852 [ + - ][ + + ]: 462 : CubitVector ctrl_pts[6];
1853 : : int jj;
1854 [ + - ][ + + ]: 396 : for ( ii = myEdgeList.size(); ii > 0; ii-- )
1855 : : {
1856 [ + - ]: 330 : edge_ptr = myEdgeList.get_and_step();
1857 [ + - ]: 330 : edge_ptr->control_points(ctrl_pts);
1858 [ + + ]: 1320 : for (jj=1; jj<4; jj++)
1859 : : {
1860 [ + - ]: 990 : x = ctrl_pts[jj].x();
1861 [ + - ]: 990 : y = ctrl_pts[jj].y();
1862 [ + - ]: 990 : z = ctrl_pts[jj].z();
1863 [ - + ]: 990 : if ( x < x_min )
1864 : 0 : x_min = x;
1865 [ - + ]: 990 : if ( y < y_min )
1866 : 0 : y_min = y;
1867 [ - + ]: 990 : if ( z < z_min )
1868 : 0 : z_min = z;
1869 [ - + ]: 990 : if ( x > x_max )
1870 : 0 : x_max = x;
1871 [ - + ]: 990 : if ( y > y_max )
1872 : 0 : y_max = y;
1873 [ - + ]: 990 : if ( z > z_max )
1874 : 0 : z_max = z;
1875 : : }
1876 : : }
1877 : : CubitFacet *facet_ptr;
1878 [ + - ][ + + ]: 198 : for ( ii = myFacetList.size(); ii > 0; ii-- )
1879 : : {
1880 [ + - ]: 132 : facet_ptr = myFacetList.get_and_step();
1881 [ + - ]: 132 : facet_ptr->get_control_points(ctrl_pts);
1882 [ + + ]: 924 : for (jj=0; jj<6; jj++)
1883 : : {
1884 [ + - ]: 792 : x = ctrl_pts[jj].x();
1885 [ + - ]: 792 : y = ctrl_pts[jj].y();
1886 [ + - ]: 792 : z = ctrl_pts[jj].z();
1887 [ - + ]: 792 : if ( x < x_min )
1888 : 0 : x_min = x;
1889 [ - + ]: 792 : if ( y < y_min )
1890 : 0 : y_min = y;
1891 [ - + ]: 792 : if ( z < z_min )
1892 : 0 : z_min = z;
1893 [ - + ]: 792 : if ( x > x_max )
1894 : 0 : x_max = x;
1895 [ - + ]: 792 : if ( y > y_max )
1896 : 0 : y_max = y;
1897 [ - + ]: 792 : if ( z > z_max )
1898 : 0 : z_max = z;
1899 : : }
1900 : : }
1901 : : }
1902 [ + - ]: 605 : CubitVector min_v(x_min, y_min, z_min );
1903 [ + - ]: 605 : CubitVector max_v(x_max, y_max, z_max );
1904 [ + - ][ + - ]: 605 : myBBox = new CubitBox( min_v, max_v );
1905 : : }
1906 : 2838 : return *(myBBox);
1907 : : }
1908 : :
1909 : : //===========================================================================
1910 : : //Function Name: closest_point
1911 : : //
1912 : : //Member Type: PUBLIC
1913 : : //Description: Finds the closest point from the vector (this_point) to the
1914 : : // set of facets that lies on the set of facets. If the point
1915 : : // lies outside this set, the closest point will be on the plane
1916 : : // of the closest facet. The closest_point is set to be that point.
1917 : : //===========================================================================
1918 : 2167 : CubitStatus FacetEvalTool::closest_point(CubitVector &this_point,
1919 : : CubitVector *closest_point_ptr,
1920 : : CubitVector *normal_ptr)
1921 : : {
1922 : 2167 : CubitBoolean trim = CUBIT_FALSE;
1923 : : CubitBoolean outside;
1924 : 2167 : CubitStatus rv = CUBIT_SUCCESS;
1925 : 2167 : static int count = 0; count++;
1926 : :
1927 : 2167 : int mydebug = 0;
1928 [ - + ]: 2167 : if (mydebug)
1929 : : {
1930 [ # # ]: 0 : debug_draw_vec( this_point, CUBIT_RED_INDEX );
1931 : : }
1932 : :
1933 : : rv = project_to_facets( this_point, trim, &outside,
1934 [ + - ]: 2167 : closest_point_ptr, normal_ptr );
1935 : :
1936 [ + - ][ + - ]: 2167 : if (DEBUG_FLAG(49))
[ - + ]
1937 : : {
1938 [ # # ]: 0 : if (closest_point_ptr)
1939 : : {
1940 [ # # ]: 0 : double dist = closest_point_ptr->distance_between( this_point );
1941 [ # # ]: 0 : if (dist > 1.0)
1942 : : {
1943 [ # # ][ # # ]: 0 : PRINT_ERROR("Appears to be bad projection in FacetEvalTool::project_to_facets\n");
[ # # ][ # # ]
1944 : : }
1945 : : }
1946 : : }
1947 : :
1948 [ - + ][ # # ]: 2167 : if (mydebug && closest_point_ptr)
1949 : : {
1950 [ # # ]: 0 : debug_draw_vec( *closest_point_ptr, CUBIT_GREEN_INDEX );
1951 : : }
1952 : 2167 : return rv;
1953 : : }
1954 : :
1955 : :
1956 : 0 : CubitFacet* FacetEvalTool::closest_facet( const CubitVector& point )
1957 : : {
1958 [ # # ]: 0 : CubitVector non_const(point);
1959 : : CubitBoolean junk;
1960 [ # # ]: 0 : CubitStatus result = project_to_facets( non_const, true, &junk, 0, 0 );
1961 [ # # ]: 0 : return result ? lastFacet : 0;
1962 : : }
1963 : :
1964 : : //===========================================================================
1965 : : //Function Name: facets_from_search_grid
1966 : : //
1967 : : //Member Type: PRIVATE
1968 : : //Description: find the closest facets to the point in the search grid
1969 : : // by starting with a default tolerance and expanding until facets are found
1970 : : //===========================================================================
1971 : 0 : void FacetEvalTool::facets_from_search_grid(
1972 : : CubitVector &this_point,
1973 : : DLIList<CubitFacet *> &facet_list,
1974 : : double &tol_used)
1975 : : {
1976 : 0 : double tol = compareTol * 10;
1977 [ # # ]: 0 : while (facet_list.size() == 0)
1978 : : {
1979 [ # # ]: 0 : CubitVector ptmin( this_point.x() - tol,
1980 [ # # ]: 0 : this_point.y() - tol,
1981 [ # # ][ # # ]: 0 : this_point.z() - tol );
1982 [ # # ]: 0 : CubitVector ptmax( this_point.x() + tol,
1983 [ # # ]: 0 : this_point.y() + tol,
1984 [ # # ][ # # ]: 0 : this_point.z() + tol );
1985 [ # # ]: 0 : CubitBox ptbox(ptmin,ptmax);
1986 [ # # ]: 0 : aTree->find(ptbox,facet_list);
1987 : :
1988 [ # # ][ # # ]: 0 : if (0 == facet_list.size())
1989 : 0 : tol *= 2.0;
1990 [ # # ]: 0 : }
1991 : :
1992 : 0 : tol_used = tol;
1993 : 0 : }
1994 : :
1995 : : //===========================================================================
1996 : : //Function Name: facets_from_search_grid
1997 : : //
1998 : : //Member Type: PRIVATE
1999 : : //Description: find the closest facets to the point in the search grid
2000 : : // for a given tolerance
2001 : : //===========================================================================
2002 : 0 : void FacetEvalTool::facets_from_search_grid(
2003 : : CubitVector &this_point,
2004 : : double compare_tol,
2005 : : DLIList<CubitFacet *> &facet_list )
2006 : : {
2007 : 0 : double tol = compare_tol;
2008 : :
2009 [ # # ]: 0 : CubitVector ptmin( this_point.x() - tol,
2010 [ # # ]: 0 : this_point.y() - tol,
2011 [ # # ][ # # ]: 0 : this_point.z() - tol );
2012 [ # # ]: 0 : CubitVector ptmax( this_point.x() + tol,
2013 [ # # ]: 0 : this_point.y() + tol,
2014 [ # # ][ # # ]: 0 : this_point.z() + tol );
2015 [ # # ]: 0 : CubitBox ptbox(ptmin,ptmax);
2016 [ # # ][ # # ]: 0 : aTree->find(ptbox,facet_list);
2017 : 0 : }
2018 : :
2019 : : //===========================================================================
2020 : : //Function Name: project_to_facets
2021 : : //
2022 : : //Member Type: PRIVATE
2023 : : //Description: Project a point to the facets. Use the interpOrder.
2024 : : // if trim is set, then trim the point to the boundary.
2025 : : // This is a non-static version. it uses the facets
2026 : : // in the evaltool
2027 : : //===========================================================================
2028 : : CubitStatus
2029 : 2178 : FacetEvalTool::project_to_facets(CubitVector &this_point,
2030 : : CubitBoolean trim,
2031 : : CubitBoolean *outside,
2032 : : CubitVector *closest_point_ptr,
2033 : : CubitVector *normal_ptr)
2034 : : {
2035 [ + - ]: 2178 : CpuTimer function_time;
2036 [ + - ][ + - ]: 2178 : if (DEBUG_FLAG(110))
[ - + ]
2037 : : {
2038 [ # # ]: 0 : function_time.cpu_secs();
2039 : 0 : numEvals++;
2040 : : }
2041 : :
2042 : 2178 : CubitStatus rv = CUBIT_SUCCESS;
2043 : :
2044 : : // if there are a lot of facets on this surface - use the grid search first
2045 : : // to narrow the selection
2046 : :
2047 [ - + ]: 2178 : if (aTree != NULL)
2048 : : {
2049 [ # # ]: 0 : DLIList<CubitFacet *> facet_list;
2050 : 0 : double search_tol = DBL_MAX;
2051 [ # # ]: 0 : facets_from_search_grid( this_point, facet_list, search_tol );
2052 [ # # ][ # # ]: 0 : if ( DEBUG_FLAG(110) )
[ # # ]
2053 [ # # ]: 0 : timeGridSearch += function_time.cpu_secs();
2054 : :
2055 [ # # ][ # # ]: 0 : if (facet_list.size())
2056 : : {
2057 [ # # ]: 0 : CubitVector grid_close_pt;
2058 : : rv = project_to_facets(facet_list,lastFacet,interpOrder,compareTol,
2059 : : this_point,trim,outside, &grid_close_pt,
2060 [ # # ]: 0 : normal_ptr);
2061 : :
2062 [ # # ]: 0 : if (CUBIT_SUCCESS == rv)
2063 : : {
2064 [ # # ]: 0 : if (closest_point_ptr)
2065 : : {
2066 [ # # ]: 0 : *closest_point_ptr = grid_close_pt;
2067 : : }
2068 : :
2069 : : // when we do the projection, if we end up with the closest point being farther
2070 : : // away than the grid search tolerance, we may have missed a closer facet
2071 : : // so redo the grid search using the distance as a tolerance
2072 [ # # ]: 0 : double distance = grid_close_pt.distance_between( this_point );
2073 [ # # ]: 0 : if (distance > search_tol)
2074 : : {
2075 [ # # ]: 0 : DLIList<CubitFacet*> facets_within_distance;
2076 [ # # ]: 0 : CubitVector grid_close_pt2;
2077 [ # # ]: 0 : facets_from_search_grid(this_point, distance, facets_within_distance);
2078 [ # # ][ # # ]: 0 : if (facets_within_distance.size() &&
[ # # ][ # # ]
2079 [ # # ]: 0 : (facets_within_distance != facet_list) )
2080 : : {
2081 : : rv = project_to_facets(facets_within_distance, lastFacet, interpOrder, compareTol,
2082 : : this_point, trim, outside, &grid_close_pt2,
2083 [ # # ]: 0 : normal_ptr);
2084 [ # # ]: 0 : if (CUBIT_SUCCESS == rv)
2085 : : {
2086 [ # # ]: 0 : double distance2 = grid_close_pt2.distance_between( this_point );
2087 [ # # ]: 0 : if (distance2 < distance)
2088 : : {
2089 [ # # ]: 0 : if (closest_point_ptr)
2090 : : {
2091 [ # # ]: 0 : *closest_point_ptr = grid_close_pt2;
2092 : : }
2093 : : }
2094 : : }
2095 [ # # ]: 0 : }
2096 : : }
2097 : : }
2098 : :
2099 [ # # ][ # # ]: 0 : if ( DEBUG_FLAG(110) )
[ # # ]
2100 : : {
2101 [ # # ]: 0 : timeFacetProject += function_time.cpu_secs();
2102 [ # # ]: 0 : if (closest_point_ptr)
2103 : : {
2104 [ # # ]: 0 : double dist = closest_point_ptr->distance_between( this_point );
2105 [ # # ]: 0 : if (dist > compareTol * 100)
2106 : : {
2107 [ # # ][ # # ]: 0 : PRINT_ERROR("Appears to be bad projection in FacetEvalTool::project_to_facets\n");
[ # # ][ # # ]
2108 [ # # ]: 0 : dcolor(CUBIT_GREEN_INDEX);
2109 [ # # ]: 0 : dpoint(this_point);
2110 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
2111 [ # # ]: 0 : dpoint(*closest_point_ptr);
2112 [ # # ]: 0 : dcolor(CUBIT_YELLOW_INDEX);
2113 [ # # ]: 0 : dfldraw(facet_list);
2114 [ # # ]: 0 : dview();
2115 : 0 : rv = CUBIT_FAILURE;
2116 : : }
2117 : : }
2118 : : }
2119 [ # # ]: 0 : }
2120 : : }
2121 : :
2122 : : // otherwise just use the complete list of facets
2123 : :
2124 : : else
2125 : : {
2126 : : rv = project_to_facets(myFacetList,lastFacet,interpOrder,compareTol,
2127 : : this_point,trim,outside,closest_point_ptr,
2128 [ + - ]: 2178 : normal_ptr);
2129 [ + - ][ + - ]: 2178 : if ( DEBUG_FLAG(110) )
[ - + ]
2130 [ # # ]: 0 : timeFacetProject += function_time.cpu_secs();
2131 : : }
2132 : :
2133 : 2178 : return rv;
2134 : : }
2135 : :
2136 : : //===========================================================================
2137 : : //Function Name: project_to_facets
2138 : : //
2139 : : //Member Type: PRIVATE
2140 : : //Description: Project a point to the facets. Use the interpOrder.
2141 : : // if trim is set, then trim the point to the boundary.
2142 : : // This is a static version of the above. Any list of facets
2143 : : // can be passed to this function
2144 : : //===========================================================================
2145 : : CubitStatus
2146 : 2178 : FacetEvalTool::project_to_facets(
2147 : : DLIList <CubitFacet *> &facet_list, // (IN) facets that we can project to
2148 : : CubitFacet *&last_facet, // (IN/OUT) last facet projected to -
2149 : : // it will try this one first
2150 : : int interp_order, // (IN) 0 = linear facets,
2151 : : // 4 = b-spline patches
2152 : : double compare_tol, // (IN) tolerance for projection -
2153 : : // should be about 1e-3*edge
2154 : : const CubitVector &this_point, // (IN) point we are projecting
2155 : : CubitBoolean trim, // (IN) trim to facet (always trimmed
2156 : : // if b-spline patch)
2157 : : CubitBoolean *outside, // (OUT) TRUE if projected outside
2158 : : // the facet
2159 : : CubitVector *closest_point_ptr, // (OUT) resulting projection point
2160 : : // (NULL if only want normal)
2161 : : CubitVector *normal_ptr) // (OUT) resulting normal at projected
2162 : : // point (NULL if not required)
2163 : : {
2164 : 2178 : int ncheck, ii, nincr=0;
2165 : : static int calls=0;
2166 : : static int nncheck=0;
2167 : : static int ntol=0;
2168 : : static int mydebug=0;
2169 : : CubitBoolean outside_facet, best_outside_facet;
2170 [ + - ][ + - ]: 2178 : CubitVector close_point, best_point, best_areacoord;
[ + - ]
2171 : 2178 : CubitFacet *best_facet = NULL;
2172 : : CubitFacet *facet;
2173 [ + - ][ - + ]: 2178 : assert (facet_list.size() > 0);
2174 : 2178 : double big_dist = compare_tol * 1.0e3;
2175 : :
2176 : : // set the first facet to be checked as the last one located
2177 : :
2178 [ + + ]: 2178 : if (last_facet) {
2179 [ + - ][ - + ]: 1958 : if (!facet_list.move_to(last_facet)) {
2180 [ # # ]: 0 : facet_list.reset();
2181 : : }
2182 : : }
2183 : : else {
2184 [ + - ]: 220 : facet_list.reset();
2185 : : }
2186 : :
2187 : : // so we don't evaluate a facet more than once - mark the facets
2188 : : // as we evaluate them. Put the evaluated facets on a used_facet_list
2189 : : // so we clear the marks off when we are done. Note: this assumes
2190 : : // theat marks are initially cleared.
2191 : :
2192 [ + - ]: 2178 : DLIList<CubitFacet *>used_facet_list;
2193 [ + - ][ + + ]: 10890 : for(ii=0; ii<facet_list.size(); ii++)
2194 [ + - ][ + - ]: 8712 : facet_list.get_and_step()->marked(0);
2195 : :
2196 [ + - ]: 2178 : int nfacets = facet_list.size();
2197 : 2178 : int nevald = 0;
2198 : 2178 : double tol = compare_tol * 10;
2199 : 2178 : const double atol = 0.001;
2200 : 2178 : double mindist = CUBIT_DBL_MAX;
2201 : 2178 : CubitBoolean eval_all = CUBIT_FALSE;
2202 : 2178 : CubitBoolean done = CUBIT_FALSE;
2203 : 2178 : best_outside_facet = CUBIT_TRUE;
2204 : :
2205 [ + + ]: 7062 : while(!done) {
2206 : :
2207 : : // define a bounding box around the point
2208 : :
2209 [ + - ]: 4884 : double ptmin_x = this_point.x() - tol;
2210 [ + - ]: 4884 : double ptmin_y = this_point.y() - tol;
2211 [ + - ]: 4884 : double ptmin_z = this_point.z() - tol;
2212 [ + - ]: 4884 : double ptmax_x = this_point.x() + tol;
2213 [ + - ]: 4884 : double ptmax_y = this_point.y() + tol;
2214 [ + - ]: 4884 : double ptmax_z = this_point.z() + tol;
2215 : :
2216 : 4884 : bool debug = false;
2217 : :
2218 [ - + ]: 4884 : if( debug )
2219 : : {
2220 [ # # ]: 0 : GfxDebug::clear();
2221 [ # # ][ # # ]: 0 : for( int i=used_facet_list.size(); i--; )
2222 [ # # ][ # # ]: 0 : used_facet_list.get_and_step()->debug_draw(CUBIT_GREEN_INDEX, 0 );
2223 : : }
2224 : :
2225 : 4884 : ncheck = 0;
2226 [ + - ][ + + ]: 19833 : for ( ii = facet_list.size(); ii > 0 && !done; ii-- )
[ + + ]
2227 : : {
2228 [ + - ]: 14949 : facet = facet_list.get_and_step();
2229 [ + - ][ + + ]: 14949 : if (facet->marked())
2230 : 11121 : continue;
2231 : :
2232 : : // Try to trivially reject this facet with a bounding box test
2233 : : // (Does the bounding box of the facet intersect with the
2234 : : // bounding box of the point)
2235 : :
2236 [ - + ]: 14553 : if( debug )
2237 : : {
2238 [ # # ]: 0 : facet->debug_draw( CUBIT_RED_INDEX );
2239 [ # # ]: 0 : GfxDebug::flush();
2240 [ # # ]: 0 : GfxDebug::mouse_xforms();
2241 : : }
2242 : :
2243 : :
2244 [ + - ]: 14553 : if (!eval_all)
2245 : : {
2246 [ + - ]: 14553 : const CubitBox &bbox = facet->bounding_box();
2247 [ + - ][ + + ]: 23111 : if (ptmax_x < bbox.min_x() ||
[ + + ][ + + ]
2248 [ + - ]: 8558 : ptmin_x > bbox.max_x()) {
2249 : 7227 : continue;
2250 : : }
2251 [ + - ][ + + ]: 12826 : if (ptmax_y < bbox.min_y() ||
[ + + ][ + + ]
2252 [ + - ]: 5500 : ptmin_y > bbox.max_y()) {
2253 : 2607 : continue;
2254 : : }
2255 [ + - ][ + + ]: 8822 : if (ptmax_z < bbox.min_z() ||
[ + + ][ + + ]
2256 [ + - ]: 4103 : ptmin_z > bbox.max_z()) {
2257 : 4719 : continue;
2258 : : }
2259 : : }
2260 : :
2261 : : // skip zero area facets
2262 [ + - ][ - + ]: 3828 : if(facet->area() <= 0.0)
2263 : : {
2264 [ # # ]: 0 : facet->marked(1);
2265 : 0 : nfacets--;
2266 : 0 : continue;
2267 : : }
2268 : :
2269 : : // Only facets that pass the bounding box test will get past here!
2270 : :
2271 : : // Project point to plane of the facet and determine its area coordinates
2272 : :
2273 : 3828 : ncheck++; nevald++;
2274 [ + - ]: 3828 : CubitVector pt_on_plane;
2275 : : double dist_to_plane;
2276 [ + - ]: 3828 : project_to_facet_plane( facet, this_point, pt_on_plane, dist_to_plane );
2277 : :
2278 [ + - ]: 3828 : CubitVector areacoord;
2279 [ + - ]: 3828 : facet_area_coordinate( facet, pt_on_plane, areacoord );
2280 [ - + ]: 3828 : if (interp_order != 0)
2281 : : {
2282 : :
2283 : : // modify the areacoord - project to the bezier patch- snaps to the
2284 : : // edge of the patch if necessary
2285 : :
2286 [ # # ]: 0 : if (project_to_facet( facet, this_point, areacoord,
2287 [ # # ]: 0 : close_point, outside_facet, compare_tol ) != CUBIT_SUCCESS)
2288 : : {
2289 : 0 : return CUBIT_FAILURE;
2290 : : }
2291 : : }
2292 : : else
2293 : : {
2294 : :
2295 : : // If sign of areacoords are all positive then its inside the triangle
2296 : : // and we are done - go interpolate the point. (use an absolute
2297 : : // tolerance since the coordinates arenormalized)
2298 [ + - ][ + + ]: 10516 : if (areacoord.x() > -atol &&
[ + + ]
2299 [ + + ][ + - ]: 6688 : areacoord.y() > -atol &&
[ + + ]
2300 [ + - ]: 2211 : areacoord.z() > -atol) {
2301 [ + + ][ + + ]: 2079 : if (dist_to_plane < compare_tol && !trim)
2302 : : {
2303 [ + - ]: 1738 : close_point = this_point;
2304 : : }
2305 : : else
2306 : : {
2307 [ + - ]: 341 : close_point = pt_on_plane;
2308 : : }
2309 : 2079 : outside_facet = CUBIT_FALSE;
2310 : : }
2311 : :
2312 : : // otherwise find the closest vertex or edge to the projected point
2313 : :
2314 [ + - ][ + + ]: 1749 : else if (areacoord.x() < atol)
2315 : : {
2316 : 968 : outside_facet = CUBIT_TRUE;
2317 [ + - ][ + + ]: 968 : if (areacoord.y() < atol)
2318 : : {
2319 [ + - ][ - + ]: 429 : if (eval_point( facet, 2, close_point ) != CUBIT_SUCCESS)
2320 : : {
2321 : 0 : return CUBIT_FAILURE;
2322 : : }
2323 : : }
2324 [ + - ][ + + ]: 539 : else if(areacoord.z() < atol)
2325 : : {
2326 [ + - ][ - + ]: 44 : if (eval_point( facet, 1, close_point ) != CUBIT_SUCCESS)
2327 : : {
2328 : 0 : return CUBIT_FAILURE;
2329 : : }
2330 : : }
2331 : : else
2332 : : {
2333 [ - + ]: 495 : if(project_to_facetedge( facet, 1, 2, this_point, pt_on_plane,
2334 [ + - ]: 495 : close_point, outside_facet, trim ) !=CUBIT_SUCCESS)
2335 : : {
2336 : 0 : return CUBIT_FAILURE;
2337 : : }
2338 : : }
2339 : : }
2340 [ + - ][ + + ]: 781 : else if (areacoord.y() < atol)
2341 : : {
2342 : 649 : outside_facet = CUBIT_TRUE;
2343 [ + - ][ + + ]: 649 : if (areacoord.z() < atol)
2344 : : {
2345 [ + - ][ - + ]: 44 : if (eval_point( facet, 0, close_point )!= CUBIT_SUCCESS)
2346 : : {
2347 : 0 : return CUBIT_FAILURE;
2348 : : }
2349 : : }
2350 : : else
2351 : : {
2352 [ - + ]: 605 : if(project_to_facetedge( facet, 2, 0, this_point, pt_on_plane,
2353 [ + - ]: 605 : close_point, outside_facet, trim ) !=CUBIT_SUCCESS)
2354 : : {
2355 : 0 : return CUBIT_FAILURE;
2356 : : }
2357 : : }
2358 : : }
2359 : : else
2360 : : {
2361 : 132 : outside_facet = CUBIT_TRUE;
2362 [ - + ]: 132 : if(project_to_facetedge( facet, 0, 1, this_point, pt_on_plane,
2363 [ + - ]: 132 : close_point, outside_facet, trim ) !=CUBIT_SUCCESS)
2364 : : {
2365 : 0 : return CUBIT_FAILURE;
2366 : : }
2367 : : }
2368 : : }
2369 : :
2370 : : // keep track of the minimum distance
2371 [ + - ]: 3828 : double dist = close_point.distance_between( this_point );
2372 : :
2373 [ + + ]: 3828 : if( trim )
2374 : : {
2375 [ + + ]: 44 : if( dist < mindist )
2376 : : {
2377 : 11 : mindist = dist;
2378 [ + - ]: 11 : best_point = close_point;
2379 : 11 : best_facet = facet;
2380 [ + - ]: 11 : best_areacoord = areacoord;
2381 : 11 : best_outside_facet = outside_facet;
2382 : : }
2383 : : }
2384 [ + + ][ + + ]: 3784 : else if ( (best_outside_facet == outside_facet && dist < mindist) ||
[ + + ]
2385 [ + + ][ - + ]: 2706 : (best_outside_facet && !outside_facet && (dist < big_dist || !best_facet)) )
[ # # ]
2386 : : {
2387 : 2662 : mindist = dist;
2388 [ + - ]: 2662 : best_point = close_point;
2389 : 2662 : best_facet = facet;
2390 [ + - ]: 2662 : best_areacoord = areacoord;
2391 : 2662 : best_outside_facet = outside_facet;
2392 : :
2393 [ + + ][ + - ]: 2662 : if (dist < compare_tol && !trim) {
2394 : 1738 : done = CUBIT_TRUE;
2395 : : }
2396 : 2662 : big_dist = 10.0 * mindist;
2397 : : }
2398 [ + - ]: 3828 : facet->marked(1);
2399 [ + - ]: 3828 : used_facet_list.append(facet);
2400 : : }
2401 : :
2402 : : // We are done if we found at least one triangle. Otherwise
2403 : : // increase the tolerance and try again
2404 : :
2405 [ + + ]: 4884 : if (!done)
2406 : : {
2407 [ + + ]: 3146 : if (nevald == nfacets)
2408 : : {
2409 : 440 : done = CUBIT_TRUE;
2410 : : }
2411 : : else
2412 : : {
2413 : 2706 : nincr++;
2414 [ + + ]: 2706 : if (ncheck > 0) {
2415 [ + - ]: 132 : if (best_outside_facet) {
2416 [ + - ]: 132 : if (nincr < 10)
2417 : : {
2418 : 132 : tol *= 2.0;
2419 : 132 : ntol++;
2420 : : }
2421 : : else
2422 : : // getting here means that the compare_tol probably is too small
2423 : : // just try all the remaining facets
2424 : : {
2425 : 132 : eval_all = CUBIT_TRUE;
2426 : : }
2427 : : }
2428 : : else
2429 : : {
2430 : 132 : done = CUBIT_TRUE;
2431 : : }
2432 : : }
2433 : : else {
2434 [ - + ]: 2574 : if (nincr >= 10)
2435 : : {
2436 : 0 : eval_all = CUBIT_TRUE;
2437 : : }
2438 : : else
2439 : : {
2440 : 2574 : tol *= 2.0e0;
2441 : 3146 : ntol++;
2442 : : }
2443 : : }
2444 : : }
2445 : : }
2446 : : } // while(!done)
2447 [ - + ]: 2178 : if(best_facet == NULL){
2448 [ # # ][ # # ]: 0 : PRINT_ERROR("Unable to determine facet correctly.\n");
[ # # ][ # # ]
2449 : 0 : return CUBIT_FAILURE;
2450 : : }
2451 : : // make sure we actually got something
2452 [ - + ]: 2178 : assert(best_facet != NULL);
2453 : :
2454 : : // if the closest point is outside of a facet, then evaluate the point
2455 : : // on the facet using its area coordinates (otherwise it would be
2456 : : // trimmed to an edge or point)
2457 : :
2458 : :
2459 [ + + ][ + + ]: 2178 : if ( !trim && best_outside_facet && interp_order != 4) {
[ + - ]
2460 [ - + ]: 132 : if (project_to_facet( best_facet, this_point, best_areacoord,
2461 [ + - ]: 132 : best_point, best_outside_facet, compare_tol )
2462 : : != CUBIT_SUCCESS) {
2463 : 0 : return CUBIT_FAILURE;
2464 : : }
2465 : :
2466 : : // see if its really outside (it could just be on an edge where the
2467 : : // curvature is convex)
2468 : :
2469 [ + - ]: 132 : best_outside_facet = is_outside( best_facet, best_areacoord );
2470 : : }
2471 : :
2472 : : // evaluate the normal if required
2473 : :
2474 [ + + ]: 2178 : if (normal_ptr) {
2475 [ + - ]: 924 : CubitVector normal;
2476 [ + - ][ - + ]: 924 : if (eval_facet_normal( best_facet, best_areacoord, normal )
2477 : : != CUBIT_SUCCESS) {
2478 : 0 : return CUBIT_FAILURE;
2479 : : }
2480 [ + - ]: 924 : *normal_ptr = normal;
2481 : : }
2482 : :
2483 [ + + ]: 2178 : if (closest_point_ptr) {
2484 [ + - ]: 1276 : *closest_point_ptr = best_point;
2485 : : }
2486 : :
2487 : 2178 : *outside = best_outside_facet;
2488 : 2178 : last_facet = best_facet;
2489 : :
2490 : :
2491 : : // clear the marks from the used facets
2492 : :
2493 [ + - ][ + + ]: 6006 : for (ii=0; ii<used_facet_list.size(); ii++)
2494 : : {
2495 [ + - ]: 3828 : facet = used_facet_list.get_and_step();
2496 [ + - ]: 3828 : facet->marked( 0 );
2497 : : }
2498 : :
2499 [ - + ]: 2178 : if (mydebug) {
2500 : 0 : nncheck+= ncheck;
2501 : 0 : calls++;
2502 [ # # ]: 0 : if (calls%100==0){
2503 [ # # ][ # # ]: 0 : PRINT_INFO("calls = %d, ckecks = %d, ntol = %d\n",calls,nncheck,ntol);
[ # # ][ # # ]
2504 : : }
2505 : : }
2506 : :
2507 [ + - ]: 2178 : return CUBIT_SUCCESS;
2508 : : }
2509 : :
2510 : :
2511 : : //===========================================================================
2512 : : //Function Name: eval_facet
2513 : : //
2514 : : //Member Type: PUBLIC
2515 : : //Descriptoin: Evaluate the location and normal of a set of area coordinates
2516 : : // on a facet. Use the interpOrder to evaluate.
2517 : : // Static function
2518 : : //===========================================================================
2519 : 0 : CubitStatus FacetEvalTool::eval_facet( CubitFacet *facet,
2520 : : CubitVector &areacoord,
2521 : : CubitVector *eval_point,
2522 : : CubitVector *eval_normal )
2523 : : {
2524 : 0 : CubitStatus rv = CUBIT_SUCCESS;
2525 : 0 : int mydebug = 0;
2526 [ # # ]: 0 : if (mydebug)
2527 : : {
2528 [ # # ]: 0 : dcolor( CUBIT_RED_INDEX );
2529 [ # # ]: 0 : dfdraw( facet );
2530 [ # # ]: 0 : dview();
2531 [ # # ]: 0 : dcolor( CUBIT_YELLOW_INDEX );
2532 : : }
2533 : :
2534 [ # # ]: 0 : CubitPoint *pt0 = facet->point(0);
2535 [ # # ]: 0 : CubitPoint *pt1 = facet->point(1);
2536 [ # # ]: 0 : CubitPoint *pt2 = facet->point(2);
2537 [ # # ]: 0 : CubitVector close_point;
2538 : :
2539 [ # # ]: 0 : CubitVector this_point;
2540 [ # # ][ # # ]: 0 : this_point.x( areacoord.x() * pt0->x() +
2541 [ # # ][ # # ]: 0 : areacoord.y() * pt1->x() +
2542 [ # # ][ # # ]: 0 : areacoord.z() * pt2->x() );
[ # # ]
2543 [ # # ][ # # ]: 0 : this_point.y( areacoord.x() * pt0->y() +
2544 [ # # ][ # # ]: 0 : areacoord.y() * pt1->y() +
2545 [ # # ][ # # ]: 0 : areacoord.z() * pt2->y() );
[ # # ]
2546 [ # # ][ # # ]: 0 : this_point.z( areacoord.x() * pt0->z() +
2547 [ # # ][ # # ]: 0 : areacoord.y() * pt1->z() +
2548 [ # # ][ # # ]: 0 : areacoord.z() * pt2->z() );
[ # # ]
2549 : :
2550 [ # # ]: 0 : int eval_order = facet->eval_order();
2551 [ # # ]: 0 : if (eval_order != 0)
2552 : : {
2553 [ # # ][ # # ]: 0 : if (facet->is_flat())
2554 : 0 : eval_order = 0;
2555 : : }
2556 : :
2557 [ # # # ]: 0 : switch(eval_order) {
2558 : : case 0:
2559 [ # # ]: 0 : if (eval_point)
2560 [ # # ]: 0 : *eval_point = this_point;
2561 [ # # ]: 0 : if (eval_normal)
2562 [ # # ]: 0 : eval_facet_normal(facet, areacoord, *eval_normal);
2563 : 0 : break;
2564 : :
2565 : : case 4:
2566 : : eval_bezier_patch( facet,
2567 : : areacoord,
2568 [ # # ]: 0 : close_point );
2569 : : //project_to_patch(facet, areacoord, this_point, close_point,
2570 : : // eval_normal, outside );
2571 : :
2572 : : //for now over-ride the normal from project_to_patch -- it is bogus
2573 [ # # ]: 0 : if (eval_normal)
2574 [ # # ]: 0 : eval_facet_normal(facet, areacoord, *eval_normal);
2575 [ # # ]: 0 : if (eval_point)
2576 [ # # ]: 0 : *eval_point = close_point;
2577 : :
2578 : 0 : break;
2579 : : default:
2580 : : // the interpolation order for now is limited to 0 and 4
2581 : : // something other than that is being attempted (or eval_order is not
2582 : : // returning the correct value)
2583 : 0 : assert(0);
2584 : : break;
2585 : : }
2586 : :
2587 : 0 : return rv;
2588 : : }
2589 : :
2590 : : //===========================================================================
2591 : : //Function Name: eval_facet
2592 : : //
2593 : : //Member Type: PUBLIC
2594 : : //Descriptoin: Evaluate the location of a set of area coordinates
2595 : : // on a facet. Use the interpOrder to evaluate.
2596 : : //===========================================================================
2597 : 0 : CubitStatus FacetEvalTool::eval_facet( CubitFacet *facet,
2598 : : CubitVector &pt,
2599 : : CubitVector &areacoord,
2600 : : CubitVector &close_point,
2601 : : CubitBoolean &outside_facet )
2602 : : {
2603 : :
2604 : 0 : CubitPoint *pt0 = facet->point(0);
2605 : 0 : CubitPoint *pt1 = facet->point(1);
2606 : 0 : CubitPoint *pt2 = facet->point(2);
2607 : :
2608 : 0 : int interp_order = facet->eval_order();
2609 [ # # ][ # # ]: 0 : if (interp_order != 0 && facet->is_flat())
[ # # ]
2610 : : {
2611 : 0 : interp_order = 0;
2612 : : }
2613 : :
2614 [ # # # # : 0 : switch(interp_order) {
# # ]
2615 : : case 0:
2616 : 0 : close_point.x( areacoord.x() * pt0->x() +
2617 : 0 : areacoord.y() * pt1->x() +
2618 : 0 : areacoord.z() * pt2->x() );
2619 : 0 : close_point.y( areacoord.x() * pt0->y() +
2620 : 0 : areacoord.y() * pt1->y() +
2621 : 0 : areacoord.z() * pt2->y() );
2622 : 0 : close_point.z( areacoord.x() * pt0->z() +
2623 : 0 : areacoord.y() * pt1->z() +
2624 : 0 : areacoord.z() * pt2->z() );
2625 : 0 : outside_facet = CUBIT_FALSE;
2626 : 0 : break;
2627 : : case 1:
2628 : : {
2629 [ # # ]: 0 : CubitVector tp0 = pt0->project_to_tangent_plane( pt );
2630 [ # # ]: 0 : CubitVector tp1 = pt1->project_to_tangent_plane( pt );
2631 [ # # ]: 0 : CubitVector tp2 = pt2->project_to_tangent_plane( pt );
2632 [ # # ][ # # ]: 0 : close_point.x( areacoord.x() * tp0.x() +
2633 [ # # ][ # # ]: 0 : areacoord.y() * tp1.x() +
2634 [ # # ][ # # ]: 0 : areacoord.z() * tp2.x() );
[ # # ]
2635 [ # # ][ # # ]: 0 : close_point.y( areacoord.x() * tp0.y() +
2636 [ # # ][ # # ]: 0 : areacoord.y() * tp1.y() +
2637 [ # # ][ # # ]: 0 : areacoord.z() * tp2.y() );
[ # # ]
2638 [ # # ][ # # ]: 0 : close_point.z( areacoord.x() * tp0.z() +
2639 [ # # ][ # # ]: 0 : areacoord.y() * tp1.z() +
2640 [ # # ][ # # ]: 0 : areacoord.z() * tp2.z() );
[ # # ]
2641 : 0 : outside_facet = CUBIT_FALSE;
2642 : : }
2643 : 0 : break;
2644 : : case 2:
2645 : : {
2646 [ # # ][ # # ]: 0 : CubitVector qp0, qp1, qp2, qn0, qn1, qn2;
[ # # ][ # # ]
[ # # ][ # # ]
2647 [ # # ]: 0 : eval_quadratic( facet, 0, pt, qp0, qn0 );
2648 [ # # ]: 0 : eval_quadratic( facet, 1, pt, qp1, qn1 );
2649 [ # # ]: 0 : eval_quadratic( facet, 2, pt, qp2, qn2 );
2650 : :
2651 [ # # ][ # # ]: 0 : close_point.x( areacoord.x() * qp0.x() +
2652 [ # # ][ # # ]: 0 : areacoord.y() * qp1.x() +
2653 [ # # ][ # # ]: 0 : areacoord.z() * qp2.x() );
[ # # ]
2654 [ # # ][ # # ]: 0 : close_point.y( areacoord.x() * qp0.y() +
2655 [ # # ][ # # ]: 0 : areacoord.y() * qp1.y() +
2656 [ # # ][ # # ]: 0 : areacoord.z() * qp2.y() );
[ # # ]
2657 [ # # ][ # # ]: 0 : close_point.z( areacoord.x() * qp0.z() +
2658 [ # # ][ # # ]: 0 : areacoord.y() * qp1.z() +
2659 [ # # ][ # # ]: 0 : areacoord.z() * qp2.z() );
[ # # ]
2660 : 0 : outside_facet = CUBIT_FALSE;
2661 : : }
2662 : 0 : break;
2663 : : case 3:
2664 : : {
2665 [ # # ][ # # ]: 0 : CubitVector qp0, qp1, qp2;
[ # # ]
2666 [ # # ]: 0 : eval_quadric( facet, 0, pt, qp0 );
2667 [ # # ]: 0 : eval_quadric( facet, 1, pt, qp1 );
2668 [ # # ]: 0 : eval_quadric( facet, 2, pt, qp2 );
2669 : :
2670 [ # # ][ # # ]: 0 : close_point.x( areacoord.x() * qp0.x() +
2671 [ # # ][ # # ]: 0 : areacoord.y() * qp1.x() +
2672 [ # # ][ # # ]: 0 : areacoord.z() * qp2.x() );
[ # # ]
2673 [ # # ][ # # ]: 0 : close_point.y( areacoord.x() * qp0.y() +
2674 [ # # ][ # # ]: 0 : areacoord.y() * qp1.y() +
2675 [ # # ][ # # ]: 0 : areacoord.z() * qp2.y() );
[ # # ]
2676 [ # # ][ # # ]: 0 : close_point.z( areacoord.x() * qp0.z() +
2677 [ # # ][ # # ]: 0 : areacoord.y() * qp1.z() +
2678 [ # # ][ # # ]: 0 : areacoord.z() * qp2.z() );
[ # # ]
2679 : 0 : outside_facet = CUBIT_FALSE;
2680 : : }
2681 : 0 : break;
2682 : : case 4:
2683 : : {
2684 : : //CubitStatus stat = eval_bezier_patch( facet, areacoord, pt );
2685 : : //close_point = pt;
2686 : : //outside_facet = CUBIT_FALSE;
2687 : 0 : double compare_tol = /*sqrt(facet->area()) **/ 1.0e-3;
2688 : 0 : int edge_id = -1;
2689 : : project_to_patch(facet, areacoord, pt, close_point, NULL,
2690 : 0 : outside_facet, compare_tol, edge_id );
2691 : : }
2692 : 0 : break;
2693 : : }
2694 : :
2695 : 0 : return CUBIT_SUCCESS;
2696 : : }
2697 : :
2698 : : //===========================================================================
2699 : : //Function Name: project_to_facet
2700 : : //
2701 : : //Member Type: PUBLIC
2702 : : //Description: project to a single facet. Uses the input areacoord as
2703 : : // a starting guess.
2704 : : //===========================================================================
2705 : 132 : CubitStatus FacetEvalTool::project_to_facet( CubitFacet *facet,
2706 : : const CubitVector &pt,
2707 : : CubitVector &areacoord,
2708 : : CubitVector &close_point,
2709 : : CubitBoolean &outside_facet,
2710 : : double compare_tol)
2711 : : {
2712 : :
2713 : 132 : CubitStatus stat = CUBIT_SUCCESS;
2714 : 132 : CubitPoint *pt0 = facet->point(0);
2715 : 132 : CubitPoint *pt1 = facet->point(1);
2716 : 132 : CubitPoint *pt2 = facet->point(2);
2717 : :
2718 : 132 : int interp_order = facet->eval_order();
2719 [ - + ]: 132 : if (facet->is_flat())
2720 : : {
2721 : 0 : interp_order = 0;
2722 : : }
2723 : :
2724 [ + - - - ]: 132 : switch(interp_order) {
2725 : : case 0:
2726 : 132 : close_point.x( areacoord.x() * pt0->x() +
2727 : 132 : areacoord.y() * pt1->x() +
2728 : 132 : areacoord.z() * pt2->x() );
2729 : 132 : close_point.y( areacoord.x() * pt0->y() +
2730 : 132 : areacoord.y() * pt1->y() +
2731 : 132 : areacoord.z() * pt2->y() );
2732 : 132 : close_point.z( areacoord.x() * pt0->z() +
2733 : 132 : areacoord.y() * pt1->z() +
2734 : 132 : areacoord.z() * pt2->z() );
2735 : 132 : outside_facet = CUBIT_FALSE;
2736 : 132 : break;
2737 : : case 1:
2738 : : case 2:
2739 : : case 3:
2740 : 0 : assert(0); // not available from this function
2741 : : break;
2742 : : case 4:
2743 : : {
2744 : 0 : int edge_id = -1;
2745 : : stat = project_to_patch(facet, areacoord, pt, close_point, NULL,
2746 : 0 : outside_facet, compare_tol, edge_id );
2747 : : }
2748 : 0 : break;
2749 : : }
2750 : :
2751 : 132 : return stat;
2752 : : }
2753 : :
2754 : : //===========================================================================
2755 : : //Function Name: project_to_facetedge
2756 : : //
2757 : : //Member Type: PUBLIC
2758 : : //Description: Project a point to the facet edge.
2759 : : // Use the interpOrder to evaluate.
2760 : : //===========================================================================
2761 : 1232 : CubitStatus FacetEvalTool::project_to_facetedge( CubitFacet *facet,
2762 : : int vert0, int vert1,
2763 : : const CubitVector &the_point,
2764 : : CubitVector &pt_on_plane,
2765 : : CubitVector &close_point,
2766 : : CubitBoolean &outside_facet,
2767 : : CubitBoolean must_be_on_edge)
2768 : : {
2769 [ + - ]: 1232 : CubitPoint *pt0 = facet->point(vert0);
2770 [ + - ]: 1232 : CubitPoint *pt1 = facet->point(vert1);
2771 : :
2772 : : // the edge vector
2773 : :
2774 [ + - ][ + - ]: 1232 : CubitVector e0 ( pt1->x() - pt0->x(),
2775 [ + - ][ + - ]: 1232 : pt1->y() - pt0->y(),
2776 [ + - ][ + - ]: 2464 : pt1->z() - pt0->z() );
[ + - ]
2777 [ + - ]: 1232 : e0.normalize();
2778 : :
2779 : : // vector from vert0 to point
2780 : :
2781 [ + - ][ + - ]: 1232 : CubitVector v0 ( pt_on_plane.x() - pt0->x(),
2782 [ + - ][ + - ]: 1232 : pt_on_plane.y() - pt0->y(),
2783 [ + - ][ + - ]: 2464 : pt_on_plane.z() - pt0->z() );
[ + - ]
2784 : :
2785 [ + - ][ + - ]: 1232 : CubitVector v1 ( pt_on_plane.x() - pt1->x(),
2786 [ + - ][ + - ]: 1232 : pt_on_plane.y() - pt1->y(),
2787 [ + - ][ + - ]: 2464 : pt_on_plane.z() - pt1->z() );
[ + - ]
2788 : :
2789 : : // project to edge
2790 : :
2791 [ + - ]: 1232 : double projection1 = v0%e0;
2792 [ + - ][ + - ]: 1232 : double projection2 = v1%(-e0);
2793 : :
2794 [ - + ][ # # ]: 1232 : if( !must_be_on_edge || (projection1 > 0 && projection2 > 0 ))
[ # # ]
2795 : : {
2796 [ + - ][ + - ]: 1232 : close_point.x ( pt0->x() + e0.x() * projection1 );
[ + - ]
2797 [ + - ][ + - ]: 1232 : close_point.y ( pt0->y() + e0.y() * projection1 );
[ + - ]
2798 [ + - ][ + - ]: 1232 : close_point.z ( pt0->z() + e0.z() * projection1 );
[ + - ]
2799 : : }
2800 : : else //we are closer to one a facet vertex than to the edge
2801 : : {
2802 [ # # ][ # # ]: 0 : if( the_point.distance_between( pt0->coordinates() )
[ # # ]
2803 [ # # ][ # # ]: 0 : < the_point.distance_between( pt1->coordinates()))
2804 [ # # ][ # # ]: 0 : close_point = pt0->coordinates();
2805 : : else
2806 [ # # ][ # # ]: 0 : close_point = pt1->coordinates();
2807 : : }
2808 : :
2809 : : // project the point on the facet (if the order is higher than 0)
2810 : :
2811 : 1232 : outside_facet = CUBIT_TRUE;
2812 [ + - ][ - + ]: 1232 : if (facet->eval_order() > 0 && !facet->is_flat()) {
[ # # ][ # # ]
[ - + ]
2813 [ # # ]: 0 : CubitVector areacoord;
2814 [ # # ]: 0 : facet_area_coordinate( facet, close_point, areacoord );
2815 : 0 : int edge_id = -1;
2816 [ # # ][ # # ]: 0 : if ((vert0 == 0 && vert1 == 1) ||
[ # # ]
2817 [ # # ]: 0 : (vert0 == 1 && vert1 == 0))
2818 : : {
2819 : 0 : edge_id = 2;
2820 : : }
2821 [ # # ][ # # ]: 0 : else if ((vert0 == 1 && vert1 == 2) ||
[ # # ]
2822 [ # # ]: 0 : (vert0 == 2 && vert1 == 1))
2823 : : {
2824 : 0 : edge_id = 0;
2825 : : }
2826 [ # # ][ # # ]: 0 : else if ((vert0 == 2 && vert1 == 0) ||
[ # # ]
2827 [ # # ]: 0 : (vert0 == 0 && vert1 == 2))
2828 : : {
2829 : 0 : edge_id = 1;
2830 : : }
2831 : : else
2832 : : {
2833 : 0 : assert(0); //edge_id wasn't set
2834 : : }
2835 : :
2836 : 0 : double compare_tol = projection1 * 1.0e-3;
2837 [ # # ]: 0 : if (project_to_patch( facet, areacoord, the_point, close_point,
2838 [ # # ]: 0 : NULL, outside_facet, compare_tol, edge_id )!= CUBIT_SUCCESS) {
2839 : 0 : return CUBIT_FAILURE;
2840 : : }
2841 : : }
2842 : :
2843 : 1232 : return CUBIT_SUCCESS;
2844 : : }
2845 : :
2846 : : //===========================================================================
2847 : : //Function Name: eval_edge
2848 : : //
2849 : : //Member Type: PUBLIC
2850 : : //Description: Evaluate the location of a point projected to a
2851 : : // linear edge.
2852 : : //===========================================================================
2853 : 0 : CubitStatus FacetEvalTool::eval_edge( CubitFacet *facet,
2854 : : int vert0, int vert1,
2855 : : CubitVector &pt_on_plane,
2856 : : CubitVector &close_point )
2857 : : {
2858 [ # # ]: 0 : CubitPoint *pt0 = facet->point(vert0);
2859 [ # # ]: 0 : CubitPoint *pt1 = facet->point(vert1);
2860 : :
2861 : : // the edge vector
2862 : :
2863 [ # # ][ # # ]: 0 : CubitVector e0 ( pt1->x() - pt0->x(),
2864 [ # # ][ # # ]: 0 : pt1->y() - pt0->y(),
2865 [ # # ][ # # ]: 0 : pt1->z() - pt0->z() );
[ # # ]
2866 [ # # ]: 0 : e0.normalize();
2867 : :
2868 : : // vector from vert0 to point
2869 : :
2870 [ # # ][ # # ]: 0 : CubitVector v0 ( pt_on_plane.x() - pt0->x(),
2871 [ # # ][ # # ]: 0 : pt_on_plane.y() - pt0->y(),
2872 [ # # ][ # # ]: 0 : pt_on_plane.z() - pt0->z() );
[ # # ]
2873 : :
2874 : : // project to edge
2875 : :
2876 [ # # ]: 0 : double len = v0%e0;
2877 [ # # ][ # # ]: 0 : close_point.x ( pt0->x() + e0.x() * len );
[ # # ]
2878 [ # # ][ # # ]: 0 : close_point.y ( pt0->y() + e0.y() * len );
[ # # ]
2879 [ # # ][ # # ]: 0 : close_point.z ( pt0->z() + e0.z() * len );
[ # # ]
2880 : :
2881 : 0 : return CUBIT_SUCCESS;
2882 : : }
2883 : :
2884 : : //===========================================================================
2885 : : //Function Name: eval_point
2886 : : //
2887 : : //Member Type: PUBLIC
2888 : : //Descriptoin: Evaluate the location and normal of a set of area coordinates
2889 : : // on a facet.
2890 : : //===========================================================================
2891 : 517 : CubitStatus FacetEvalTool::eval_point( CubitFacet *facet,
2892 : : int vertex_id,
2893 : : CubitVector &close_point )
2894 : : {
2895 : 517 : close_point.x (facet->point(vertex_id)->x());
2896 : 517 : close_point.y (facet->point(vertex_id)->y());
2897 : 517 : close_point.z (facet->point(vertex_id)->z());
2898 : :
2899 : 517 : return CUBIT_SUCCESS;
2900 : : }
2901 : :
2902 : : //===========================================================================
2903 : : //Function Name: eval_facet_normal
2904 : : //
2905 : : //Member Type: PUBLIC
2906 : : //Descriptoin: Evaluate the normal of the facet (use the interpOrder)
2907 : : // return normalized normal
2908 : : //===========================================================================
2909 : 924 : CubitStatus FacetEvalTool::eval_facet_normal( CubitFacet *facet,
2910 : : CubitVector &areacoord,
2911 : : CubitVector &normal )
2912 : : {
2913 [ + - - - ]: 924 : switch(facet->eval_order()) {
2914 : : case 0:
2915 [ + - ]: 924 : normal = facet->normal();
2916 : 924 : break;
2917 : : case 1: case 2: case 3:
2918 : : {
2919 [ # # ][ # # ]: 0 : CubitVector norm0 = facet->point(0)->normal( facet );
2920 [ # # ][ # # ]: 0 : CubitVector norm1 = facet->point(1)->normal( facet );
2921 [ # # ][ # # ]: 0 : CubitVector norm2 = facet->point(2)->normal( facet );
2922 [ # # ][ # # ]: 0 : normal.x( areacoord.x() * norm0.x() +
2923 [ # # ][ # # ]: 0 : areacoord.y() * norm1.x() +
2924 [ # # ][ # # ]: 0 : areacoord.z() * norm2.x() );
[ # # ]
2925 [ # # ][ # # ]: 0 : normal.y( areacoord.x() * norm0.y() +
2926 [ # # ][ # # ]: 0 : areacoord.y() * norm1.y() +
2927 [ # # ][ # # ]: 0 : areacoord.z() * norm2.y() );
[ # # ]
2928 [ # # ][ # # ]: 0 : normal.z( areacoord.x() * norm0.z() +
2929 [ # # ][ # # ]: 0 : areacoord.y() * norm1.z() +
2930 [ # # ][ # # ]: 0 : areacoord.z() * norm2.z() );
[ # # ]
2931 [ # # ]: 0 : normal.normalize();
2932 : : }
2933 : 0 : break;
2934 : : case 4:
2935 [ # # ]: 0 : if(facet->is_flat())
2936 : : {
2937 [ # # ]: 0 : normal = facet->normal();
2938 : : }
2939 : : else
2940 : : {
2941 : 0 : eval_bezier_patch_normal(facet, areacoord, normal);
2942 : :
2943 : : // check for reasonableness of the normal
2944 : :
2945 : : #if 0
2946 : : if (DEBUG_FLAG(110))
2947 : : {
2948 : : CubitVector norm0 = facet->point(0)->normal( facet );
2949 : : CubitVector norm1 = facet->point(1)->normal( facet );
2950 : : CubitVector norm2 = facet->point(2)->normal( facet );
2951 : : CubitVector lin_normal;
2952 : : lin_normal.x( areacoord.x() * norm0.x() +
2953 : : areacoord.y() * norm1.x() +
2954 : : areacoord.z() * norm2.x() );
2955 : : lin_normal.y( areacoord.x() * norm0.y() +
2956 : : areacoord.y() * norm1.y() +
2957 : : areacoord.z() * norm2.y() );
2958 : : lin_normal.z( areacoord.x() * norm0.z() +
2959 : : areacoord.y() * norm1.z() +
2960 : : areacoord.z() * norm2.z() );
2961 : : lin_normal.normalize();
2962 : :
2963 : : PRINT_INFO("(facet %4d) ac=%7.5lf %7.5lf %7.5lf\n",
2964 : : facet->id(), areacoord.x(), areacoord.y(), areacoord.z());
2965 : : PRINT_INFO(" bn=%7.5lf %7.5lf %7.5lf\n",
2966 : : normal.x(), normal.y(), normal.z());
2967 : : PRINT_INFO(" ln=%7.5lf %7.5lf %7.5lf\n",
2968 : : lin_normal.x(), lin_normal.y(), lin_normal.z());
2969 : :
2970 : : const double tol = 1e-2;
2971 : : if (fabs(lin_normal.x() - normal.x()) > tol ||
2972 : : fabs(lin_normal.y() - normal.y()) > tol ||
2973 : : fabs(lin_normal.z() - normal.z()) > tol)
2974 : : {
2975 : : int mydebug = 0;
2976 : : if (mydebug)
2977 : : {
2978 : : CubitVector pt;
2979 : : eval_bezier_patch(facet, areacoord, pt);
2980 : : dcolor(CUBIT_GREEN_INDEX);
2981 : : dray(pt,normal,1.0);
2982 : : dcolor(CUBIT_RED_INDEX);
2983 : : dray(pt,lin_normal,1.0);
2984 : : dview();
2985 : : }
2986 : :
2987 : : PRINT_INFO("^=============^==============^=============^=============^\n");
2988 : : }
2989 : : }
2990 : : #endif
2991 : : }
2992 : :
2993 : 0 : break;
2994 : : }
2995 : :
2996 : 924 : return CUBIT_SUCCESS;
2997 : : }
2998 : :
2999 : :
3000 : : //===========================================================================
3001 : : //Function Name: eval_quadratic
3002 : : //
3003 : : //Member Type: PRIVATE
3004 : : //Descriptoin: Evaluate the point based on a quadratic approximation
3005 : : //===========================================================================
3006 : 0 : CubitStatus FacetEvalTool::eval_quadratic( CubitFacet *facet,
3007 : : int pt_idx,
3008 : : CubitVector &eval_pt,
3009 : : CubitVector &qpoint,
3010 : : CubitVector &qnorm )
3011 : : {
3012 : : // interpolate a point on a circle that is defined by two points and
3013 : : // two normals. The first pont is a vertex on the facet, the second
3014 : : // is a point on the opposite edge. The point to be interpolated lies
3015 : : // somewhere between the two
3016 : :
3017 [ # # ][ # # ]: 0 : CubitVector normA = facet->point(pt_idx)->normal( facet );
3018 [ # # ][ # # ]: 0 : CubitVector ptA = facet->point(pt_idx)->coordinates();
3019 : 0 : int idx0 = -1, idx1 = -1;
3020 [ # # # # ]: 0 : switch(pt_idx) {
3021 : : case 0:
3022 : 0 : idx0 = 1;
3023 : 0 : idx1 = 2;
3024 : 0 : break;
3025 : : case 1:
3026 : 0 : idx0 = 2;
3027 : 0 : idx1 = 0;
3028 : 0 : break;
3029 : : case 2:
3030 : 0 : idx0 = 0;
3031 : 0 : idx1 = 1;
3032 : 0 : break;
3033 : : }
3034 [ # # ][ # # ]: 0 : CubitVector ptB, normB, pt0, pt1, norm0, norm1;
[ # # ][ # # ]
[ # # ][ # # ]
3035 [ # # ][ # # ]: 0 : pt0 = facet->point(idx0)->coordinates();
[ # # ]
3036 [ # # ][ # # ]: 0 : pt1 = facet->point(idx1)->coordinates();
[ # # ]
3037 [ # # ][ # # ]: 0 : norm0 = facet->point(idx0)->normal( facet );
[ # # ]
3038 [ # # ][ # # ]: 0 : norm1 = facet->point(idx1)->normal( facet );
[ # # ]
3039 : : on_circle( pt0, norm0, pt1, norm1,
3040 [ # # ]: 0 : eval_pt, ptB, normB );
3041 : : on_circle( ptA, normA, ptB, normB,
3042 [ # # ]: 0 : eval_pt, qpoint, qnorm );
3043 : :
3044 : 0 : return CUBIT_SUCCESS;
3045 : :
3046 : : }
3047 : :
3048 : : //===========================================================================
3049 : : //Function Name: on_circle
3050 : : //
3051 : : //Member Type: PRIVATE
3052 : : //Description: compute a projection of a point onto a circle. The circle
3053 : : // is defined by two points and two normals. Return the
3054 : : // projected point and its normal
3055 : : //===========================================================================
3056 : 0 : void FacetEvalTool::on_circle( CubitVector &ptA,
3057 : : CubitVector &normA,
3058 : : CubitVector &ptB,
3059 : : CubitVector &normB,
3060 : : CubitVector &eval_pt,
3061 : : CubitVector &pt_on_circle,
3062 : : CubitVector &normal )
3063 : : {
3064 : : // angle between the normals
3065 : :
3066 : 0 : double cosang = normA%normB;
3067 : :
3068 : : // check for flat surfaces - project to the segment
3069 : :
3070 [ # # ][ # # ]: 0 : if (cosang >= 0.99999 || cosang <= -0.99999) {
3071 [ # # ]: 0 : CubitVector vAB = ptB - ptA;
3072 [ # # ]: 0 : vAB.normalize();
3073 [ # # ]: 0 : CubitVector vAeval_pt = eval_pt - ptA;
3074 [ # # ]: 0 : double len = vAB%vAeval_pt;
3075 [ # # ][ # # ]: 0 : pt_on_circle = ptA + len * vAB;
[ # # ]
3076 [ # # ]: 0 : if (cosang <= -0.99999) { // this is bad! (facet spans 180 degrees)
3077 [ # # ]: 0 : normal = normA;
3078 : : }
3079 : : else {
3080 [ # # ][ # # ]: 0 : normal = 0.5e0 * (normA + normB);
[ # # ]
3081 : 0 : }
3082 : : }
3083 : : else {
3084 : :
3085 : : // curved surface
3086 : : // define a common plane at eval_pt
3087 : :
3088 [ # # ]: 0 : CubitVector pnorm = normA*normB;
3089 [ # # ]: 0 : pnorm.normalize();
3090 [ # # ][ # # ]: 0 : double pcoefd = -pnorm%eval_pt;
3091 : :
3092 : : // project everything to common plane
3093 : :
3094 [ # # ]: 0 : double pdist = pnorm%ptA + pcoefd;
3095 [ # # ][ # # ]: 0 : CubitVector pptA = ptA - pnorm * pdist;
3096 [ # # ]: 0 : pdist = pnorm%ptB + pcoefd;
3097 [ # # ][ # # ]: 0 : CubitVector pptB = ptB - pnorm * pdist;
3098 : :
3099 : 0 : double angle = acos(cosang);
3100 [ # # ]: 0 : double dist = pptA.distance_between(pptB);
3101 : :
3102 : : // kradius is the radius of curvature
3103 : : // center is the center of curvature
3104 : : // centerA and centerB should be the same within tol
3105 : :
3106 : 0 : double kradius = dist / (2.0e0 * sin( angle * 0.5e0 ));
3107 [ # # ][ # # ]: 0 : CubitVector centerA = pptA - kradius * normA;
3108 [ # # ][ # # ]: 0 : CubitVector centerB = pptB - kradius * normB;
3109 [ # # ][ # # ]: 0 : CubitVector center = (centerA + centerB) * 0.5e0;
3110 : :
3111 [ # # ][ # # ]: 0 : normal = eval_pt - center;
3112 [ # # ]: 0 : normal.normalize();
3113 [ # # ][ # # ]: 0 : pt_on_circle = center + kradius * normal;
[ # # ]
3114 : : }
3115 : 0 : }
3116 : :
3117 : : //===========================================================================
3118 : : //Function Name: eval_quadric
3119 : : //
3120 : : //Member Type: PRIVATE
3121 : : //Description: evaluate a point on an interpolated quadric surface
3122 : : //===========================================================================
3123 : 0 : void FacetEvalTool::eval_quadric( CubitFacet *facet,
3124 : : int pt_index,
3125 : : CubitVector &eval_pt,
3126 : : CubitVector &qpt )
3127 : : {
3128 : : // transform the point to the local system
3129 : :
3130 [ # # ]: 0 : CubitPoint *point = facet->point(pt_index);
3131 [ # # ]: 0 : CubitVector loc_eval_pt;
3132 [ # # ]: 0 : point->transform_to_local( eval_pt, loc_eval_pt );
3133 : : // point->transform_to_local( point->coordinates(), loc_pt );
3134 : :
3135 : : // interpolate a "z" value in the local coordinate system
3136 : :
3137 [ # # ]: 0 : double *coef = point->coef_vector();
3138 [ # # ]: 0 : loc_eval_pt.z( coef[0] * loc_eval_pt.x() +
3139 [ # # ]: 0 : coef[1] * loc_eval_pt.y() +
3140 [ # # ][ # # ]: 0 : coef[2] * FacetEvalToolUtils::sqr(loc_eval_pt.x()) +
3141 [ # # ][ # # ]: 0 : coef[3] * loc_eval_pt.x() * loc_eval_pt.y() +
3142 [ # # ][ # # ]: 0 : coef[4] * FacetEvalToolUtils::sqr(loc_eval_pt.y()) );
[ # # ]
3143 : :
3144 : :
3145 : : // loc_eval_pt.z( point->z() -
3146 : : // coef[0] * loc_eval_pt.x() -
3147 : : // coef[1] * loc_eval_pt.y() +
3148 : : // coef[2] * sqr(loc_eval_pt.x()) +
3149 : : // coef[3] * loc_eval_pt.x() * loc_eval_pt.y() +
3150 : : // coef[4] * sqr(loc_eval_pt.y()) );
3151 : :
3152 : : // transform back to global system
3153 : :
3154 [ # # ]: 0 : point->transform_to_global( loc_eval_pt, qpt );
3155 : :
3156 : 0 : }
3157 : :
3158 : : //===========================================================================
3159 : : //Function Name: project_to_facet_plane
3160 : : //
3161 : : //Member Type: PUBLIC
3162 : : //Descriptoin: Project a point to the plane of a facet
3163 : : //===========================================================================
3164 : 3828 : void FacetEvalTool::project_to_facet_plane(
3165 : : CubitFacet *facet,
3166 : : const CubitVector &pt,
3167 : : CubitVector &point_on_plane,
3168 : : double &dist_to_plane)
3169 : : {
3170 [ + - ]: 3828 : CubitVector normal = facet->normal();
3171 [ + - ]: 3828 : normal.normalize();
3172 [ + - ]: 3828 : CubitPoint *facPoint = facet->point(0);
3173 [ + - ][ + - ]: 3828 : double coefd = -(normal.x() * facPoint->x() +
3174 [ + - ][ + - ]: 3828 : normal.y() * facPoint->y() +
3175 [ + - ][ + - ]: 3828 : normal.z() * facPoint->z());
3176 [ + - ][ + - ]: 3828 : double dist = normal.x() * pt.x() +
3177 [ + - ][ + - ]: 3828 : normal.y() * pt.y() +
3178 [ + - ][ + - ]: 3828 : normal.z() * pt.z() + coefd;
3179 : 3828 : dist_to_plane = fabs(dist);
3180 : :
3181 [ + - ][ + - ]: 3828 : point_on_plane.set( pt.x() - normal.x() * dist,
3182 [ + - ][ + - ]: 3828 : pt.y() - normal.y() * dist,
3183 [ + - ][ + - ]: 7656 : pt.z() - normal.z() * dist );
[ + - ]
3184 : :
3185 : 3828 : }
3186 : :
3187 : : //===========================================================================
3188 : : //Function Name: facet_area_coordinate
3189 : : //
3190 : : //Member Type: PUBLIC
3191 : : //Descriptoin: Determine the area coordinates of a point on the plane
3192 : : // of a facet
3193 : : //===========================================================================
3194 : 3828 : void FacetEvalTool::facet_area_coordinate(
3195 : : CubitFacet *facet,
3196 : : const CubitVector &pt_on_plane,
3197 : : CubitVector &areacoord )
3198 : : {
3199 : : double area2;
3200 [ + - ]: 3828 : CubitVector normal = facet->normal();
3201 [ + - ]: 3828 : CubitPoint *pt0 = facet->point(0);
3202 [ + - ]: 3828 : CubitPoint *pt1 = facet->point(1);
3203 [ + - ]: 3828 : CubitPoint *pt2 = facet->point(2);
3204 : 3828 : double tol = GEOMETRY_RESABS * 1.e-5;
3205 [ + - ][ + - ]: 3828 : CubitVector v1( pt1->x() - pt0->x(),
3206 [ + - ][ + - ]: 3828 : pt1->y() - pt0->y(),
3207 [ + - ][ + - ]: 7656 : pt1->z() - pt0->z());//(*p1-*p0);
[ + - ]
3208 [ + - ][ + - ]: 3828 : CubitVector v2( pt2->x() - pt0->x(),
3209 [ + - ][ + - ]: 3828 : pt2->y() - pt0->y(),
3210 [ + - ][ + - ]: 7656 : pt2->z() - pt0->z());// = (*p2-*p0);
[ + - ]
3211 [ + - ][ + - ]: 3828 : area2 = (v1*v2).length_squared();
3212 [ - + ]: 3828 : if(area2 < 100 * tol){
3213 : 0 : tol = .01 * area2;
3214 : : }
3215 [ + - ][ + - ]: 3828 : CubitVector absnorm( fabs(normal.x()), fabs(normal.y()), fabs(normal.z()) );
[ + - ][ + - ]
3216 : :
3217 : : // project to the closest coordinate plane so we only have to do this in 2D
3218 : :
3219 [ + - ][ + - ]: 3828 : if (absnorm.x() >= absnorm.y() && absnorm.x() >= absnorm.z()) {
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
3220 : 2420 : area2 = FacetEvalToolUtils::determ3(pt0->y(), pt0->z(),
3221 : 2420 : pt1->y(), pt1->z(),
3222 [ + - ][ + - ]: 1210 : pt2->y(), pt2->z());
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
3223 [ - + ]: 1210 : if (fabs(area2) < tol) {
3224 [ # # ]: 0 : areacoord.set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
3225 : : }
3226 [ + - ][ + - ]: 1210 : else if ( pt_on_plane.within_tolerance( pt0->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3227 : : {
3228 [ + - ]: 99 : areacoord.set( 1.0, 0.0, 0.0 );
3229 : : }
3230 [ + - ][ + - ]: 1111 : else if ( pt_on_plane.within_tolerance( pt1->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3231 : : {
3232 [ + - ]: 77 : areacoord.set( 0.0, 1.0, 0.0 );
3233 : : }
3234 [ + - ][ + - ]: 1034 : else if ( pt_on_plane.within_tolerance( pt2->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3235 : : {
3236 [ + - ]: 264 : areacoord.set( 0.0, 0.0, 1.0 );
3237 : : }
3238 : : else {
3239 : : areacoord.x(
3240 : : FacetEvalToolUtils::determ3( pt_on_plane.y(), pt_on_plane.z(),
3241 [ + - ][ + - ]: 770 : pt1->y(), pt1->z(), pt2->y(), pt2->z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3242 : : areacoord.y(
3243 : 1540 : FacetEvalToolUtils::determ3( pt0->y(), pt0->z(),
3244 : : pt_on_plane.y(), pt_on_plane.z(),
3245 [ + - ][ + - ]: 770 : pt2->y(), pt2->z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3246 : : areacoord.z(
3247 : 3080 : FacetEvalToolUtils::determ3( pt0->y(), pt0->z(), pt1->y(), pt1->z(),
3248 [ + - ][ + - ]: 1210 : pt_on_plane.y(), pt_on_plane.z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3249 : : }
3250 : : }
3251 [ + - ][ + - ]: 2618 : else if(absnorm.y() >= absnorm.x() && absnorm.y() >= absnorm.z()) {
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ]
3252 : 2772 : area2 = FacetEvalToolUtils::determ3(pt0->x(), pt0->z(),
3253 : 2772 : pt1->x(), pt1->z(),
3254 [ + - ][ + - ]: 1386 : pt2->x(), pt2->z());
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
3255 [ - + ]: 1386 : if (fabs(area2) < tol) {
3256 [ # # ]: 0 : areacoord.set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
3257 : : }
3258 [ + - ][ + - ]: 1386 : else if ( pt_on_plane.within_tolerance( pt0->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3259 : : {
3260 [ + - ]: 55 : areacoord.set( 1.0, 0.0, 0.0 );
3261 : : }
3262 [ + - ][ + - ]: 1331 : else if ( pt_on_plane.within_tolerance( pt1->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3263 : : {
3264 [ + - ]: 121 : areacoord.set( 0.0, 1.0, 0.0 );
3265 : : }
3266 [ + - ][ + - ]: 1210 : else if ( pt_on_plane.within_tolerance( pt2->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3267 : : {
3268 [ + - ]: 286 : areacoord.set( 0.0, 0.0, 1.0 );
3269 : : }
3270 : : else {
3271 : : areacoord.x(
3272 : : FacetEvalToolUtils::determ3( pt_on_plane.x(), pt_on_plane.z(),
3273 [ + - ][ + - ]: 924 : pt1->x(), pt1->z(), pt2->x(), pt2->z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3274 : : areacoord.y(
3275 : 1848 : FacetEvalToolUtils::determ3( pt0->x(), pt0->z(),
3276 : : pt_on_plane.x(), pt_on_plane.z(),
3277 [ + - ][ + - ]: 924 : pt2->x(), pt2->z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3278 : : areacoord.z(
3279 : 3696 : FacetEvalToolUtils::determ3( pt0->x(), pt0->z(), pt1->x(), pt1->z(),
3280 [ + - ][ + - ]: 1386 : pt_on_plane.x(), pt_on_plane.z() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3281 : : }
3282 : : }
3283 : : else {
3284 : 2464 : area2 = FacetEvalToolUtils::determ3(pt0->x(), pt0->y(),
3285 : 2464 : pt1->x(), pt1->y(),
3286 [ + - ][ + - ]: 1232 : pt2->x(), pt2->y());
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
3287 [ - + ]: 1232 : if (fabs(area2) < tol) {
3288 [ # # ]: 0 : areacoord.set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
3289 : : }
3290 [ + - ][ + - ]: 1232 : else if ( pt_on_plane.within_tolerance( pt0->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3291 : : {
3292 [ + - ]: 110 : areacoord.set( 1.0, 0.0, 0.0 );
3293 : : }
3294 [ + - ][ + - ]: 1122 : else if ( pt_on_plane.within_tolerance( pt1->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3295 : : {
3296 [ + - ]: 66 : areacoord.set( 0.0, 1.0, 0.0 );
3297 : : }
3298 [ + - ][ + - ]: 1056 : else if ( pt_on_plane.within_tolerance( pt2->coordinates(), GEOMETRY_RESABS ) )
[ + + ]
3299 : : {
3300 [ + - ]: 440 : areacoord.set( 0.0, 0.0, 1.0 );
3301 : : }
3302 : : else {
3303 : : areacoord.x(
3304 : : FacetEvalToolUtils::determ3( pt_on_plane.x(), pt_on_plane.y(),
3305 [ + - ][ + - ]: 616 : pt1->x(), pt1->y(), pt2->x(), pt2->y() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3306 : : areacoord.y(
3307 : 1232 : FacetEvalToolUtils::determ3( pt0->x(), pt0->y(),
3308 : : pt_on_plane.x(), pt_on_plane.y(),
3309 [ + - ][ + - ]: 616 : pt2->x(), pt2->y() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3310 : : areacoord.z(
3311 : 2464 : FacetEvalToolUtils::determ3( pt0->x(), pt0->y(), pt1->x(), pt1->y(),
3312 [ + - ][ + - ]: 616 : pt_on_plane.x(), pt_on_plane.y() ) / area2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
3313 : : }
3314 : : }
3315 : 3828 : }
3316 : :
3317 : : //===========================================================================
3318 : : //Function Name: is_outside
3319 : : //
3320 : : //Member Type: PRIVATE
3321 : : //Descriptoin: Determines if the areacoord is actually outside the range
3322 : : // of the surface facets.
3323 : : //===========================================================================
3324 : 132 : CubitBoolean FacetEvalTool::is_outside(
3325 : : CubitFacet *facet,
3326 : : CubitVector &areacoord)
3327 : : {
3328 [ - + ]: 132 : if (areacoord.x() < -GEOMETRY_RESABS) {
3329 [ # # ]: 0 : if (NULL == facet->shared_facet_on_surf( facet->point(1),
3330 : 0 : facet->point(2),
3331 : 0 : facet->tool_id() )) {
3332 : 0 : return CUBIT_TRUE;
3333 : : }
3334 : : }
3335 [ + + ]: 132 : if (areacoord.y() < -GEOMETRY_RESABS) {
3336 [ - + ]: 22 : if (NULL == facet->shared_facet_on_surf( facet->point(2),
3337 : 22 : facet->point(0),
3338 : 22 : facet->tool_id() )) {
3339 : 0 : return CUBIT_TRUE;
3340 : : }
3341 : : }
3342 [ + + ]: 132 : if (areacoord.z() < -GEOMETRY_RESABS) {
3343 [ + - ]: 110 : if (NULL == facet->shared_facet_on_surf( facet->point(0),
3344 : 110 : facet->point(1),
3345 : 110 : facet->tool_id() )) {
3346 : 110 : return CUBIT_TRUE;
3347 : : }
3348 : : }
3349 : 22 : return CUBIT_FALSE;
3350 : : }
3351 : :
3352 : : //===========================================================================
3353 : : //Function Name: closest_point_trimmed
3354 : : //
3355 : : //Member Type: PUBLIC
3356 : : //Descriptoin: Finds the closest point from the vector (this_point) to the
3357 : : // set of facets that lies on the set of facets. If the point
3358 : : // lies outside this set, the closest point will be on the edge
3359 : : // of the closest facet.
3360 : : // This function also determines if the point is outside or
3361 : : // inside the current set of facets. You should call
3362 : : // a bounding box test first before this...
3363 : : //===========================================================================
3364 : 11 : CubitStatus FacetEvalTool::closest_point_trimmed(
3365 : : CubitVector &this_point,
3366 : : CubitVector *closest_point_ptr,
3367 : : CubitBoolean &lies_outside,
3368 : : CubitVector *normal_ptr)
3369 : : {
3370 : :
3371 : 11 : CubitBoolean trim = CUBIT_TRUE;
3372 : : return project_to_facets ( this_point, trim, &lies_outside,
3373 : 11 : closest_point_ptr, normal_ptr );
3374 : : }
3375 : :
3376 : : //===========================================================================
3377 : : //Function Name: destroy_facets
3378 : : //
3379 : : //Member Type: PRIVATE
3380 : : //Descriptoin: Deletes the points and facets.
3381 : : //===========================================================================
3382 : 0 : void FacetEvalTool::destroy_facets()
3383 : : {
3384 : : int i, j;
3385 : : CubitPoint *point;
3386 : : CubitFacet *facet;
3387 : : CubitFacetEdge *edge;
3388 : 0 : myFacetList.last();
3389 : 0 : CubitFacetEdgeData *cfe_data = NULL;
3390 [ # # ]: 0 : for( i = myFacetList.size(); i > 0; i-- )
3391 : : {
3392 : 0 : facet = myFacetList.pop();
3393 [ # # ]: 0 : for (j = 0; j<3; j++)
3394 : : {
3395 : 0 : point = facet->point( j );
3396 : 0 : point->remove_facet( facet );
3397 : 0 : edge = facet->edge( j );
3398 [ # # ]: 0 : cfe_data = CAST_TO( edge, CubitFacetEdgeData );
3399 [ # # ]: 0 : if (cfe_data)
3400 : 0 : cfe_data->remove_facet( facet );
3401 : : }
3402 [ # # ]: 0 : delete facet;
3403 : : }
3404 [ # # ]: 0 : for( i = myEdgeList.size(); i > 0; i-- )
3405 : : {
3406 : 0 : edge = myEdgeList.pop();
3407 [ # # ][ # # ]: 0 : if (edge && edge->num_adj_facets() == 0)
[ # # ]
3408 : : {
3409 [ # # ]: 0 : delete edge;
3410 : : }
3411 : : }
3412 [ # # ]: 0 : for( i = myPointList.size(); i > 0; i-- )
3413 : : {
3414 : 0 : point = myPointList.pop();
3415 [ # # ][ # # ]: 0 : if (point && point->num_adj_facets() == 0)
[ # # ]
3416 : : {
3417 [ # # ]: 0 : delete point;
3418 : : }
3419 : : }
3420 : 0 : }
3421 : :
3422 : : //===========================================================================
3423 : : //Function Name: get_points_from_facets
3424 : : //
3425 : : //Member Type: PRIVATE
3426 : : //Descriptoin: Gets the set of points contained in the list of facets.
3427 : : // Populates the point_list with those points.
3428 : : //===========================================================================
3429 : 0 : CubitStatus FacetEvalTool::get_points_from_facets(
3430 : : DLIList<CubitFacet*> &facet_list,
3431 : : DLIList<CubitPoint*> &point_list )
3432 : : {
3433 : : int i, j;
3434 : :
3435 [ # # ]: 0 : DLIList<CubitPoint*> all_points;
3436 : :
3437 [ # # ][ # # ]: 0 : for ( i = 0; i < facet_list.size(); i++)
3438 : : {
3439 [ # # ]: 0 : CubitFacet* facet = facet_list.get_and_step();
3440 [ # # ]: 0 : facet->points(all_points);
3441 : : }
3442 [ # # ][ # # ]: 0 : for ( i = 0; i < all_points.size(); i++)
3443 : : {
3444 [ # # ][ # # ]: 0 : all_points.get_and_step()->marked( CUBIT_FALSE );
3445 : : }
3446 : :
3447 [ # # ][ # # ]: 0 : for ( j = 0; j < all_points.size(); j++)
3448 : : {
3449 [ # # ]: 0 : CubitPoint* point = all_points.get_and_step();
3450 [ # # ][ # # ]: 0 : if (!point->marked())
3451 : : {
3452 [ # # ]: 0 : point->marked(CUBIT_TRUE);
3453 [ # # ]: 0 : point_list.append(point);
3454 : : }
3455 : : }
3456 : :
3457 : : // unmark the found points
3458 : :
3459 [ # # ][ # # ]: 0 : for (i = 0; i < point_list.size(); i++)
3460 : : {
3461 [ # # ][ # # ]: 0 : point_list.get_and_step()->marked(CUBIT_FALSE);
3462 : : }
3463 [ # # ]: 0 : return CUBIT_SUCCESS;
3464 : : }
3465 : :
3466 : : //===========================================================================
3467 : : //Function Name: get_loops_from_facets
3468 : : //
3469 : : //Member Type: PRIVATE
3470 : : //Descriptoin: generate an ordered list of facetedge lists representing the
3471 : : // boundary loops for this set of facets
3472 : : //===========================================================================
3473 : 341 : CubitStatus FacetEvalTool::get_loops_from_facets(
3474 : : DLIList<CubitFacetEdge*> &all_edge_list, // all the edges on the facets
3475 : : DLIList<DLIList<CubitFacetEdge*>*> &loop_list ) // return a list of edge lists
3476 : : {
3477 : : int i;
3478 : 341 : int mydebug = 0;
3479 : : CubitFacetEdge* edge, *startedge;
3480 : : CubitFacet *facet, *nextfacet;
3481 : : CubitBoolean found;
3482 : :
3483 [ - + ]: 341 : if (mydebug)
3484 : : {
3485 [ # # ]: 0 : GfxDebug::clear();
3486 [ # # ][ # # ]: 0 : for (i = 0; i< myFacetList.size(); i++)
3487 [ # # ][ # # ]: 0 : myFacetList.get_and_step()->debug_draw( CUBIT_GREEN_INDEX, false);
3488 [ # # ]: 0 : GfxDebug::flush();
3489 [ # # ]: 0 : GfxDebug::mouse_xforms();
3490 : : }
3491 : :
3492 [ + - ][ + + ]: 2563 : for ( i = 0; i < all_edge_list.size(); i++)
3493 : : {
3494 [ + - ]: 2222 : startedge = all_edge_list.get_and_step();
3495 [ + - ][ + + ]: 2222 : if (startedge->get_flag() == 0) {
3496 [ + - ]: 858 : startedge->set_flag( 1 );
3497 : :
3498 : : // Find an edge without a neighboring facet or a facet that
3499 : : // is not on the current surface to start a loop
3500 : :
3501 [ + - ][ + + ]: 858 : if (startedge->num_adj_facets_on_surf(toolID) == 1)
3502 : : {
3503 : :
3504 : : // Start a new loop
3505 : :
3506 [ + - ][ + - ]: 330 : DLIList<CubitFacetEdge*> *edge_list = new DLIList<CubitFacetEdge*>;
3507 [ + - ]: 330 : loop_list.append( edge_list );
3508 [ + - ]: 330 : edge_list->append( startedge );
3509 [ - + ][ # # ]: 330 : if (mydebug) debug_draw_edge(startedge);
3510 : :
3511 : : // find the next ccw edge on the loop. Edges are placed on the
3512 : : // list based on ccw order wrt the orientation of the facets
3513 : : // (ie. same orientation as facets)
3514 : :
3515 : 330 : edge = startedge;
3516 [ + - ]: 330 : facet = edge->adj_facet_on_surf(toolID);
3517 : 330 : int num_failures = 0;
3518 [ + + ]: 1276 : do {
3519 : 1276 : found = CUBIT_FALSE;
3520 [ + + ]: 2244 : do {
3521 [ + - ]: 2244 : edge = facet->next_edge( edge );
3522 [ - + ]: 2244 : assert(edge != NULL);
3523 [ + - ]: 2244 : edge->set_flag( 1 );
3524 : :
3525 [ + - ]: 2244 : nextfacet = edge->other_facet_on_surf(facet);
3526 : :
3527 [ + + ]: 2244 : if (nextfacet == NULL) {
3528 [ + + ]: 1276 : if (edge != startedge) {
3529 : 946 : num_failures = 0;
3530 [ + - ]: 946 : edge_list->append( edge );
3531 : : }
3532 : 1276 : found = CUBIT_TRUE;
3533 : : }
3534 : : else {
3535 : 968 : num_failures++;
3536 : 968 : facet = nextfacet;
3537 : : }
3538 [ + + ][ + - ]: 2244 : } while (!found && num_failures < myFacetList.size() );
[ + - ]
3539 [ + + ][ + - ]: 1276 : } while (edge != startedge && num_failures < myFacetList.size() );
[ + - ]
3540 : :
3541 [ - + ]: 330 : if( edge != startedge )
3542 : 858 : return CUBIT_FAILURE;
3543 : : }
3544 : : }
3545 : : }
3546 : :
3547 : : // reset the flags
3548 : :
3549 [ + - ][ + + ]: 2563 : for ( i = 0; i < all_edge_list.size(); i++) {
3550 [ + - ][ + - ]: 2222 : all_edge_list.get_and_step()->set_flag( 0 );
3551 : : }
3552 : 341 : return CUBIT_SUCCESS;
3553 : : }
3554 : :
3555 : : //===========================================================================
3556 : : //Function Name: get_edges_from_facets
3557 : : //
3558 : : //Member Type: PRIVATE
3559 : : //Descriptoin: Populates the edge list from the list of facets
3560 : : //===========================================================================
3561 : 341 : CubitStatus FacetEvalTool::get_edges_from_facets(
3562 : : DLIList<CubitFacet*> &facet_list,
3563 : : DLIList<CubitFacetEdge*> &edge_list )
3564 : : {
3565 : : int i, j, k;
3566 : : CubitPoint *p0, *p1;
3567 : : CubitFacet *facet, *adj_facet;
3568 : : CubitFacetEdge *edge;
3569 [ + - ][ + + ]: 1397 : for ( i = 0; i < facet_list.size(); i++) {
3570 [ + - ]: 1056 : facet = facet_list.get_and_step();
3571 [ + + ]: 4224 : for (j=0; j<3; j++) {
3572 [ + - ][ - + ]: 3168 : if (!(edge = facet->edge(j)))
3573 : : {
3574 [ # # ]: 0 : facet->get_edge_pts( j, p0, p1 );
3575 : 0 : edge = NULL;
3576 : 0 : k = -1;
3577 : : //find another facet on this surface sharing these two points
3578 [ # # ]: 0 : adj_facet = facet->shared_facet( p0, p1 );
3579 : :
3580 [ # # ]: 0 : if (adj_facet) {
3581 [ # # ]: 0 : edge = adj_facet->edge_from_pts(p0, p1, k);
3582 : : }
3583 : :
3584 [ # # ]: 0 : if (!edge) {
3585 : : edge = (CubitFacetEdge *)
3586 [ # # ][ # # ]: 0 : new CubitFacetEdgeData( p0, p1, facet, adj_facet, j, k );
3587 : : }
3588 : : else{
3589 [ # # ]: 0 : facet->add_edge( edge );
3590 : : }
3591 : : }
3592 [ + - ]: 3168 : edge->set_flag( 0 );
3593 : : }
3594 : : }
3595 : : // create a unique list of edges
3596 : :
3597 [ + - ][ + + ]: 1397 : for ( i = 0; i < facet_list.size(); i++)
3598 : : {
3599 [ + - ]: 1056 : facet = facet_list.get_and_step();
3600 [ + + ]: 4224 : for (j=0; j<3; j++)
3601 : : {
3602 [ + - ]: 3168 : edge = facet->edge(j);
3603 [ - + ]: 3168 : if ( !edge )
3604 : 0 : return CUBIT_FAILURE;
3605 [ + - ][ + + ]: 3168 : if (0 == edge->get_flag())
3606 : : {
3607 [ + - ]: 2222 : edge->set_flag( 1 );
3608 [ + - ]: 2222 : edge_list.append( edge );
3609 : : }
3610 : : }
3611 : : }
3612 : :
3613 : : // reset the flags on the edges
3614 : :
3615 [ + - ][ + + ]: 2563 : for ( i = 0; i < edge_list.size(); i++)
3616 : : {
3617 [ + - ]: 2222 : edge = edge_list.get_and_step();
3618 [ + - ]: 2222 : edge->set_flag( 0 );
3619 : : }
3620 : 341 : return CUBIT_SUCCESS;
3621 : : }
3622 : :
3623 : :
3624 : : //===========================================================================
3625 : : //Function Name: check_faceting
3626 : : //
3627 : : //Member Type: PRIVATE
3628 : : //Descriptoin: check the edge/face orientations
3629 : : //===========================================================================
3630 : 0 : void FacetEvalTool::check_faceting()
3631 : : {
3632 : : int ii;
3633 [ # # ][ # # ]: 0 : for (ii = 0; ii< myFacetList.size(); ii++)
3634 [ # # ][ # # ]: 0 : myFacetList.get_and_step()->debug_draw( CUBIT_YELLOW_INDEX );
3635 [ # # ]: 0 : GfxDebug::flush();
3636 [ # # ]: 0 : GfxDebug::mouse_xforms();
3637 : :
3638 [ # # ]: 0 : CubitFacet *facet = myFacetList.get_and_step();
3639 [ # # ]: 0 : CubitVector snorm = facet->normal();
3640 [ # # ]: 0 : snorm.normalize();
3641 [ # # ][ # # ]: 0 : for (ii=1; ii<myFacetList.size(); ii++)
3642 : : {
3643 [ # # ]: 0 : facet = myFacetList.get_and_step();
3644 [ # # ]: 0 : CubitVector norm = facet->normal();
3645 [ # # ]: 0 : norm.normalize();
3646 [ # # ][ # # ]: 0 : if (norm % snorm < 0.8)
3647 : : {
3648 [ # # ]: 0 : facet->debug_draw( CUBIT_RED_INDEX );
3649 : : }
3650 : : else
3651 : : {
3652 [ # # ]: 0 : facet->debug_draw( CUBIT_BLUE_INDEX );
3653 : : }
3654 : : }
3655 : 0 : }
3656 : :
3657 : :
3658 : : //===========================================================================
3659 : : //Function Name: draw_facets
3660 : : //
3661 : : //Member Type: PRIVATE
3662 : : //Descriptoin: draw the facets for debug
3663 : : //===========================================================================
3664 : 0 : void FacetEvalTool::debug_draw_facets( int color )
3665 : : {
3666 : : int ii;
3667 [ # # ]: 0 : if ( color == -1 )
3668 : 0 : color = CUBIT_YELLOW_INDEX;
3669 : 0 : CubitBoolean flush_it = CUBIT_FALSE;
3670 [ # # ]: 0 : for ( ii = myFacetList.size(); ii > 0; ii-- )
3671 : : {
3672 : 0 : myFacetList.get_and_step()->debug_draw(color, flush_it);
3673 : : }
3674 : 0 : GfxDebug::flush();
3675 : 0 : }
3676 : :
3677 : : //===========================================================================
3678 : : //Function Name: draw_vec
3679 : : //
3680 : : //Member Type: PRIVATE
3681 : : //Descriptoin: draw a single point
3682 : : //===========================================================================
3683 : 0 : void FacetEvalTool::debug_draw_vec( CubitVector &vec, int color )
3684 : : {
3685 [ # # ]: 0 : if ( color == -1 )
3686 : 0 : color = CUBIT_YELLOW_INDEX;
3687 : 0 : GfxDebug::draw_point(vec, color);
3688 : 0 : GfxDebug::flush();
3689 : 0 : }
3690 : :
3691 : : //===========================================================================
3692 : : //Function Name: draw_facet_normals
3693 : : //
3694 : : //Member Type: PRIVATE
3695 : : //Descriptoin: the the normal at the centroid of the facets
3696 : : //===========================================================================
3697 : 0 : void FacetEvalTool::debug_draw_facet_normals(int color)
3698 : : {
3699 : : int ii,jj;
3700 [ # # ]: 0 : if ( color == -1 )
3701 : 0 : color = CUBIT_RED_INDEX;
3702 [ # # ]: 0 : for ( ii = myFacetList.size(); ii > 0; ii-- )
3703 : : {
3704 [ # # ]: 0 : CubitFacet *facet = myFacetList.get_and_step();
3705 [ # # ]: 0 : CubitVector center(0.0e0, 0.0e0, 0.0e0);
3706 [ # # ]: 0 : for (jj = 0; jj < 3; jj++)
3707 : : {
3708 [ # # ][ # # ]: 0 : center += facet->point(jj)->coordinates();
[ # # ]
3709 : : }
3710 [ # # ]: 0 : center /= 3.0;
3711 [ # # ]: 0 : CubitVector normal = facet->normal();
3712 [ # # ][ # # ]: 0 : double mag = sqrt(FacetEvalToolUtils::sqr(normal.x()) + FacetEvalToolUtils::sqr(normal.y()) + FacetEvalToolUtils::sqr(normal.z()));
[ # # ][ # # ]
[ # # ][ # # ]
3713 : 0 : mag = sqrt(mag);
3714 [ # # ]: 0 : normal.normalize();
3715 [ # # ][ # # ]: 0 : CubitVector end = center + normal*mag;
3716 [ # # ][ # # ]: 0 : GfxDebug::draw_line(center.x(), center.y(), center.z(),
[ # # ]
3717 [ # # ][ # # ]: 0 : end.x(), end.y(), end.z(), color);
[ # # ][ # # ]
3718 : : }
3719 : 0 : GfxDebug::flush();
3720 : 0 : }
3721 : :
3722 : : //===========================================================================
3723 : : //Function Name: draw_point_normals
3724 : : //
3725 : : //Member Type: PRIVATE
3726 : : //Descriptoin: the the normal at the points
3727 : : //===========================================================================
3728 : 0 : void FacetEvalTool::debug_draw_point_normals(int color)
3729 : : {
3730 : : int ii;
3731 [ # # ]: 0 : if ( color == -1 )
3732 : 0 : color = CUBIT_RED_INDEX;
3733 : 0 : double len = 0.1 * sqrt(FacetEvalToolUtils::sqr(myBBox->x_range()) +
3734 : 0 : FacetEvalToolUtils::sqr(myBBox->y_range()) +
3735 : 0 : FacetEvalToolUtils::sqr(myBBox->z_range()));
3736 [ # # ]: 0 : for ( ii = myPointList.size(); ii > 0; ii-- )
3737 : : {
3738 [ # # ]: 0 : CubitPoint *point = myPointList.get_and_step();
3739 [ # # ]: 0 : CubitVector normal = point->normal();
3740 [ # # ][ # # ]: 0 : CubitVector end = point->coordinates() + normal*len;
[ # # ]
3741 [ # # ][ # # ]: 0 : GfxDebug::draw_line(point->x(), point->y(), point->z(),
[ # # ]
3742 [ # # ][ # # ]: 0 : end.x(), end.y(), end.z(), color);
[ # # ][ # # ]
3743 : : }
3744 : 0 : GfxDebug::flush();
3745 : 0 : }
3746 : :
3747 : : //===========================================================================
3748 : : //Function Name: draw_edges
3749 : : //
3750 : : //Member Type: PRIVATE
3751 : : //Descriptoin: draw the facet edges
3752 : : //===========================================================================
3753 : 0 : void FacetEvalTool::debug_draw_edges(int color)
3754 : : {
3755 : : int ii;
3756 [ # # ]: 0 : if ( color == -1 )
3757 : 0 : color = CUBIT_BLUE_INDEX;
3758 [ # # ]: 0 : for ( ii = myEdgeList.size(); ii > 0; ii-- )
3759 : : {
3760 : 0 : CubitFacetEdge *edge = myEdgeList.get_and_step();
3761 : 0 : CubitPoint *begin_point = edge->point(0);
3762 : 0 : CubitPoint *end_point = edge->point(1);
3763 : 0 : GfxDebug::draw_line(begin_point->x(),
3764 : 0 : begin_point->y(),
3765 : 0 : begin_point->z(),
3766 : 0 : end_point->x(),
3767 : 0 : end_point->y(),
3768 : 0 : end_point->z(),
3769 : 0 : color);
3770 : : }
3771 : 0 : GfxDebug::flush();
3772 : 0 : }
3773 : :
3774 : : //===========================================================================
3775 : : //Function Name: draw_edge
3776 : : //
3777 : : //Member Type: PRIVATE
3778 : : //Descriptoin: draw the facet edge
3779 : : //===========================================================================
3780 : 0 : void FacetEvalTool::debug_draw_edge(CubitFacetEdge *edge, int color)
3781 : : {
3782 [ # # ]: 0 : if ( color == -1 ) {
3783 : 0 : color = CUBIT_YELLOW_INDEX;
3784 : : }
3785 : 0 : CubitPoint *begin_point = edge->point(0);
3786 : 0 : CubitPoint *end_point = edge->point(1);
3787 : 0 : GfxDebug::draw_line(begin_point->x(),
3788 : 0 : begin_point->y(),
3789 : 0 : begin_point->z(),
3790 : 0 : end_point->x(),
3791 : 0 : end_point->y(),
3792 : 0 : end_point->z(),
3793 : 0 : color);
3794 : 0 : GfxDebug::flush();
3795 : 0 : }
3796 : :
3797 : : //===========================================================================
3798 : : //Function Name: draw_bezier_edges
3799 : : //
3800 : : //Member Type: PRIVATE
3801 : : //Descriptoin: draw the control polygons from the bezier control points on
3802 : : // the edges
3803 : : //===========================================================================
3804 : 0 : void FacetEvalTool::debug_draw_bezier_edges(int color)
3805 : : {
3806 : : int ii, i;
3807 [ # # ][ # # ]: 0 : CubitVector ctrl_pts[5], begin, end;
[ # # ][ # # ]
3808 [ # # ]: 0 : if ( color == -1 )
3809 : 0 : color = CUBIT_RED_INDEX;
3810 [ # # ][ # # ]: 0 : for ( ii = myEdgeList.size(); ii > 0; ii-- )
3811 : : {
3812 [ # # ]: 0 : CubitFacetEdge *edge = myEdgeList.get_and_step();
3813 [ # # ]: 0 : edge->control_points( ctrl_pts );
3814 [ # # ]: 0 : for (i=1; i<5; i++) {
3815 [ # # ]: 0 : begin = ctrl_pts[i-1];
3816 [ # # ]: 0 : end = ctrl_pts[i];
3817 [ # # ][ # # ]: 0 : GfxDebug::draw_line(begin.x(), begin.y(), begin.z(),
[ # # ]
3818 [ # # ][ # # ]: 0 : end.x(), end.y(), end.z(), color);
[ # # ][ # # ]
3819 : : }
3820 : : }
3821 [ # # ]: 0 : GfxDebug::flush();
3822 : 0 : }
3823 : :
3824 : : //===========================================================================
3825 : : //Function Name: draw_bezier_facets
3826 : : //
3827 : : //Member Type: PRIVATE
3828 : : //Descriptoin: draw the control polygons from the bezier control points on
3829 : : // the facets
3830 : : //===========================================================================
3831 : 0 : void FacetEvalTool::debug_draw_bezier_facets(int color)
3832 : : {
3833 : : int ii;
3834 [ # # ]: 0 : if ( color == -1 )
3835 : 0 : color = CUBIT_WHITE_INDEX;
3836 [ # # ]: 0 : for ( ii = myFacetList.size(); ii > 0; ii-- )
3837 : : {
3838 : 0 : debug_draw_bezier_facet( myFacetList.get_and_step(), color );
3839 : : }
3840 : :
3841 : 0 : }
3842 : :
3843 : : //===========================================================================
3844 : : //Function Name: draw_bezier_facet
3845 : : //
3846 : : //Member Type: PRIVATE
3847 : : //Descriptoin: draw the control polygons from the bezier control points on
3848 : : // a facet
3849 : : //===========================================================================
3850 : 0 : void FacetEvalTool::debug_draw_bezier_facet(CubitFacet *facet, int color)
3851 : : {
3852 : : CubitVector areacoord( 0.3333333333333333333,
3853 : : 0.3333333333333333333,
3854 [ # # ]: 0 : 0.3333333333333333333 );
3855 : :
3856 : : // interpolate internal control points
3857 : :
3858 [ # # ][ # # ]: 0 : CubitVector gctrl_pts[6];
3859 [ # # ]: 0 : facet->get_control_points( gctrl_pts );
3860 [ # # ][ # # ]: 0 : CubitVector P_facet[3];
3861 : :
3862 : : //2,1,1
3863 [ # # ][ # # ]: 0 : P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
[ # # ]
3864 [ # # ][ # # ]: 0 : (areacoord.y() * gctrl_pts[3] +
[ # # ]
3865 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[4]);
[ # # ]
3866 : : //1,2,1
3867 [ # # ][ # # ]: 0 : P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
[ # # ]
3868 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[0] +
[ # # ]
3869 [ # # ][ # # ]: 0 : areacoord.z() * gctrl_pts[5]);
[ # # ]
3870 : : //1,1,2
3871 [ # # ][ # # ]: 0 : P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
[ # # ]
3872 [ # # ][ # # ]: 0 : (areacoord.x() * gctrl_pts[1] +
[ # # ]
3873 [ # # ][ # # ]: 0 : areacoord.y() * gctrl_pts[2]);
[ # # ]
3874 : :
3875 [ # # ][ # # ]: 0 : CubitVector cp0[5], cp1[5], cp2[5];
[ # # ][ # # ]
[ # # ][ # # ]
3876 : :
3877 [ # # ][ # # ]: 0 : facet->edge(0)->control_points(facet,cp0);
3878 [ # # ][ # # ]: 0 : facet->edge(1)->control_points(facet,cp1);
3879 [ # # ][ # # ]: 0 : facet->edge(2)->control_points(facet,cp2);
3880 : :
3881 : 0 : color = CUBIT_WHITE_INDEX;
3882 [ # # ]: 0 : debug_draw_line(cp0[1], cp2[3], color);
3883 [ # # ]: 0 : debug_draw_line(cp0[2], P_facet[1], color);
3884 [ # # ]: 0 : debug_draw_line(P_facet[1], cp2[2], color);
3885 [ # # ]: 0 : debug_draw_line(cp0[3], P_facet[2], color);
3886 [ # # ]: 0 : debug_draw_line(P_facet[2], P_facet[0], color);
3887 [ # # ]: 0 : debug_draw_line(P_facet[0], cp2[1], color);
3888 : :
3889 : 0 : color = CUBIT_GREEN_INDEX;
3890 [ # # ]: 0 : debug_draw_line(cp1[1], P_facet[2], color);
3891 [ # # ]: 0 : debug_draw_line(P_facet[2], P_facet[1], color);
3892 [ # # ]: 0 : debug_draw_line(P_facet[1], cp2[3], color);
3893 [ # # ]: 0 : debug_draw_line(cp1[2], P_facet[0], color);
3894 [ # # ]: 0 : debug_draw_line(P_facet[0], cp2[2], color);
3895 [ # # ]: 0 : debug_draw_line(cp1[3], cp2[1], color);
3896 : :
3897 : 0 : color = CUBIT_BLUE_INDEX;
3898 [ # # ]: 0 : debug_draw_line(cp0[1], P_facet[1], color);
3899 [ # # ]: 0 : debug_draw_line(P_facet[1], P_facet[0], color);
3900 [ # # ]: 0 : debug_draw_line(P_facet[0], cp1[3], color);
3901 [ # # ]: 0 : debug_draw_line(cp0[2], P_facet[2], color);
3902 [ # # ]: 0 : debug_draw_line(P_facet[2], cp1[2], color);
3903 [ # # ]: 0 : debug_draw_line(cp0[3], cp1[1], color);
3904 : 0 : }
3905 : :
3906 : : //===========================================================================
3907 : : //Function Name: draw_eval_bezier_facet
3908 : : //
3909 : : //Member Type: PRIVATE
3910 : : //Descriptoin: draw points on the evaluated bezier patch
3911 : : //===========================================================================
3912 : 0 : void FacetEvalTool::debug_draw_eval_bezier_facet( CubitFacet *facet )
3913 : : {
3914 [ # # ][ # # ]: 0 : CubitVector areacoord, pt, loc;
[ # # ]
3915 : : double u, v, w;
3916 : : CubitBoolean outside;
3917 : : #if 0
3918 : : for (int i=0; i<=10; i++) {
3919 : : v = w = 0.5 * (double)i/10.0;
3920 : : u = 1.0 - v - w;
3921 : : areacoord.set(u, v, w);
3922 : : eval_facet(facet, pt, areacoord, loc, outside);
3923 : : draw_location(loc);
3924 : : v = 0.25 * (double)i/10.0;
3925 : : w = 0.75 * (double)i/10.0;
3926 : : u = 1.0 - v - w;
3927 : : areacoord.set(u, v, w);
3928 : : eval_facet(facet, pt, areacoord, loc, outside);
3929 : : draw_location(loc);
3930 : : w = 0.25 * (double)i/10.0;
3931 : : v = 0.75 * (double)i/10.0;
3932 : : u = 1.0 - v - w;
3933 : : areacoord.set(u, v, w);
3934 : : eval_facet(facet, pt, areacoord, loc, outside);
3935 : : debug_draw_location(loc);
3936 : : }
3937 : : #endif
3938 [ # # ]: 0 : for (int j=0; j<=10; j++) {
3939 [ # # ]: 0 : for (int i=0; i<=20; i++) {
3940 : 0 : u = ((double)j/10.0) * (double)i/20.0;
3941 : 0 : w = (1.0 - ((double)j/10.0)) * (double)i/20.0;
3942 : 0 : v = 1.0 - u - w;
3943 [ # # ]: 0 : areacoord.set(u, v, w);
3944 [ # # ]: 0 : eval_facet(facet, pt, areacoord, loc, outside);
3945 [ # # ]: 0 : debug_draw_location(loc);
3946 : : }
3947 : : }
3948 : :
3949 : : #if 0
3950 : : for (i=0; i<=10; i++) {
3951 : : u = v = 0.5 * (double)i/10.0;
3952 : : w = 1.0 - u - v;
3953 : : areacoord.set(u, v, w);
3954 : : eval_facet(facet, pt, areacoord, loc, outside);
3955 : : draw_location(loc);
3956 : : u = 0.25 * (double)i/10.0;
3957 : : v = 0.75 * (double)i/10.0;
3958 : : w = 1.0 - u - v;
3959 : : areacoord.set(u, v, w);
3960 : : eval_facet(facet, pt, areacoord, loc, outside);
3961 : : draw_location(loc);
3962 : : v = 0.25 * (double)i/10.0;
3963 : : u = 0.75 * (double)i/10.0;
3964 : : w = 1.0 - u - v;
3965 : : areacoord.set(u, v, w);
3966 : : eval_facet(facet, pt, areacoord, loc, outside);
3967 : : draw_location(loc);
3968 : : }
3969 : : #endif
3970 : :
3971 : 0 : }
3972 : :
3973 : : //===========================================================================
3974 : : //Function Name: draw_line
3975 : : //
3976 : : //Member Type: PRIVATE
3977 : : //===========================================================================
3978 : 0 : void FacetEvalTool::debug_draw_line(CubitVector &begin, CubitVector &end, int color)
3979 : : {
3980 : 0 : GfxDebug::draw_line(begin.x(), begin.y(), begin.z(),
3981 : 0 : end.x(), end.y(), end.z(), color);
3982 : 0 : GfxDebug::flush();
3983 : 0 : }
3984 : :
3985 : : //===========================================================================
3986 : : //Function Name: draw_location
3987 : : //
3988 : : //Member Type: PRIVATE
3989 : : //===========================================================================
3990 : 0 : void FacetEvalTool::debug_draw_location(CubitVector &loc, int color )
3991 : : {
3992 [ # # ]: 0 : if ( color == -1 )
3993 : 0 : color = CUBIT_YELLOW_INDEX;
3994 : 0 : GfxDebug::draw_point(loc, color);
3995 : 0 : GfxDebug::flush();
3996 : 0 : }
3997 : :
3998 : :
3999 : : //===========================================================================
4000 : : //Function Name: write_loops
4001 : : //
4002 : : //Member Type: PRIVATE
4003 : : //===========================================================================
4004 : 0 : void FacetEvalTool::write_loops()
4005 : : {
4006 : : int ii, jj;
4007 [ # # ]: 0 : for (ii=0; ii<myLoopList.size(); ii++)
4008 : : {
4009 : 0 : DLIList<CubitFacetEdge*> *loop = myLoopList.get_and_step();
4010 [ # # ][ # # ]: 0 : PRINT_INFO("======= Loop %d =========\n", ii);
4011 [ # # ]: 0 : for (jj=0; jj<loop->size(); jj++)
4012 : : {
4013 : 0 : CubitFacetEdge *edge = loop->get_and_step();
4014 : 0 : CubitPoint *point0 = edge->point( 0 );
4015 : 0 : CubitPoint *point1 = edge->point( 1 );
4016 [ # # ]: 0 : PRINT_INFO(" (%d) %f, %f, %f (%d) %f, %f, %f\n",
4017 : : point0->id(), point0->x(), point0->y(), point0->z(),
4018 [ # # ]: 0 : point1->id(), point1->x(), point1->y(), point1->z());
4019 : : }
4020 : : }
4021 : 0 : }
4022 : :
4023 : : //===========================================================================
4024 : : //Function Name: transform_control_points
4025 : : //
4026 : : //Member Type: PUBLIC
4027 : : //===========================================================================
4028 : 0 : void FacetEvalTool::transform_control_points( CubitTransformMatrix &tfmat )
4029 : : {
4030 [ # # ]: 0 : if (interpOrder != 4)
4031 : 0 : return;
4032 : :
4033 [ # # ][ # # ]: 0 : CubitVector control_points[6];
4034 : : CubitFacet *facet;
4035 : : int ii, jj;
4036 [ # # ][ # # ]: 0 : for (ii=0; ii<myFacetList.size(); ii++)
4037 : : {
4038 [ # # ]: 0 : facet = myFacetList.get_and_step();
4039 [ # # ]: 0 : facet->get_control_points( control_points );
4040 [ # # ]: 0 : for (jj=0; jj<6; jj++)
4041 : : {
4042 [ # # ][ # # ]: 0 : control_points[jj] = tfmat * control_points[jj];
4043 : : }
4044 [ # # ]: 0 : facet->set_control_points( control_points );
4045 : : }
4046 : :
4047 : : CubitFacetEdge *edge;
4048 [ # # ][ # # ]: 0 : for (ii=0; ii<myEdgeList.size(); ii++)
4049 : : {
4050 [ # # ]: 0 : edge = myEdgeList.get_and_step();
4051 [ # # ][ # # ]: 0 : if (!edge->get_flag())
4052 : : {
4053 [ # # ]: 0 : edge->set_flag(1);
4054 [ # # ]: 0 : edge->control_points( control_points );
4055 : :
4056 : : // assumes the end control points (the vertices) have already
4057 : : // been transformed
4058 : :
4059 [ # # ][ # # ]: 0 : control_points[0] = tfmat * control_points[1];
4060 [ # # ][ # # ]: 0 : control_points[1] = tfmat * control_points[2];
4061 [ # # ][ # # ]: 0 : control_points[2] = tfmat * control_points[3];
4062 [ # # ]: 0 : edge->control_points( control_points, 4 );
4063 : : }
4064 : : }
4065 : : }
4066 : :
4067 : : //===========================================================================
4068 : : //Function Name: area
4069 : : //
4070 : : //Member Type: PUBLIC
4071 : : //===========================================================================
4072 : 341 : double FacetEvalTool::area()
4073 : : {
4074 [ + - ]: 341 : if (myArea < 0.0)
4075 : 341 : calculate_area();
4076 : :
4077 : 341 : return myArea;
4078 : : }
4079 : : //===========================================================================
4080 : : //Function Name: calcualte_area
4081 : : //
4082 : : //Member Type: Public
4083 : : //===========================================================================
4084 : 407 : void FacetEvalTool::calculate_area()
4085 : : {
4086 : 407 : myArea = 0.0;
4087 : : int ii;
4088 : : CubitFacet *facet;
4089 [ + + ]: 1727 : for (ii=0; ii<myFacetList.size(); ii++)
4090 : : {
4091 : 1320 : facet = myFacetList.get_and_step();
4092 : 1320 : myArea += facet->area();
4093 : : }
4094 : 407 : }
4095 : :
4096 : : //===========================================================================
4097 : : //Function Name: parameterize
4098 : : //Description: compute the parameterization of the facetted representation
4099 : : //Member Type: PUBLIC
4100 : : //===========================================================================
4101 : 0 : CubitStatus FacetEvalTool::parameterize()
4102 : : {
4103 : :
4104 [ # # ][ # # ]: 0 : if (myLoopList.size() != 1)
4105 : : {
4106 [ # # ][ # # ]: 0 : PRINT_WARNING("Cannot parameterize surface. Multiple loops detected\n");
[ # # ][ # # ]
4107 : 0 : isParameterized = CUBIT_FALSE;
4108 : 0 : return CUBIT_FAILURE;
4109 : : }
4110 : :
4111 : : // make arrays out of the points and facets
4112 : :
4113 [ # # ][ # # ]: 0 : double *points = new double [3 * myPointList.size()];
[ # # ]
4114 [ # # ][ # # ]: 0 : int *facets = new int [3 * myFacetList.size()];
[ # # ]
4115 [ # # ][ # # ]: 0 : if (!points || !facets)
4116 : : {
4117 [ # # ][ # # ]: 0 : PRINT_ERROR("Could not define parameterization for surface (out of memory)\n");
[ # # ][ # # ]
4118 : 0 : return CUBIT_FAILURE;
4119 : : }
4120 : : int ii;
4121 : : CubitPoint *pt;
4122 [ # # ]: 0 : myPointList.reset();
4123 [ # # ][ # # ]: 0 : for (ii=0; ii<myPointList.size(); ii++)
4124 : : {
4125 [ # # ]: 0 : pt = myPointList.get_and_step();
4126 [ # # ]: 0 : points[ii*3] = pt->x();
4127 [ # # ]: 0 : points[ii*3+1] = pt->y();
4128 [ # # ]: 0 : points[ii*3+2] = pt->z();
4129 [ # # ]: 0 : pt->marked(ii);
4130 : : }
4131 : : CubitFacet *facet;
4132 : : CubitPoint *pts[3];
4133 [ # # ][ # # ]: 0 : for (ii=0; ii<myFacetList.size(); ii++)
4134 : : {
4135 [ # # ]: 0 : facet = myFacetList.get_and_step();
4136 [ # # ]: 0 : facet->points( pts[0], pts[1], pts[2] );
4137 [ # # ]: 0 : facets[ii*3] = pts[0]->marked();
4138 [ # # ]: 0 : facets[ii*3+1] = pts[1]->marked();
4139 [ # # ]: 0 : facets[ii*3+2] = pts[2]->marked();
4140 : : }
4141 : :
4142 : : // do the parameterization
4143 : :
4144 : : // comment out for now
4145 : : // Note to sjowen: this depends on FacetParamTool and facetParamTool is a ParamTool
4146 : : // (which is in geom directory). We ned a solution that breaks that dependency.
4147 : : // FacetParamTool facetparamtool( myPointList.size(), myFacetList.size(),
4148 : : // points, facets );
4149 : : if(1)//!facetparamtool.flatten())
4150 : : {
4151 [ # # ][ # # ]: 0 : PRINT_ERROR("Surface Parameterizer Failed\n");
[ # # ][ # # ]
4152 : 0 : isParameterized = CUBIT_FALSE;
4153 : : }
4154 : : else
4155 : : {
4156 : : double *sizes = NULL;
4157 : : double *uv = NULL;//facetparamtool.get_uvs_sizing( ratio, sizes );
4158 : :
4159 : : // update the points with uv values
4160 : :
4161 : : TDFacetBoundaryPoint *td_fbp;
4162 : : CubitBoolean on_internal_boundary;
4163 : : myPointList.reset();
4164 : : for (ii=0; ii<myPointList.size(); ii++)
4165 : : {
4166 : : pt = myPointList.get_and_step();
4167 : : DLIList <CubitFacet *> facet_list;
4168 : : pt->facets_on_surf( toolID, facet_list, on_internal_boundary );
4169 : : if (on_internal_boundary)
4170 : : {
4171 : : td_fbp = TDFacetBoundaryPoint::get_facet_boundary_point( pt );
4172 : : if (!td_fbp)
4173 : : {
4174 : : TDFacetBoundaryPoint::add_facet_boundary_point( pt );
4175 : : td_fbp = TDFacetBoundaryPoint::get_facet_boundary_point( pt );
4176 : : td_fbp->add_surf_facets( facet_list );
4177 : : td_fbp->set_uv( facet_list.get(), uv[ii*2], uv[ii*2+1] );
4178 : : }
4179 : : else
4180 : : {
4181 : : if (td_fbp->set_uv( facet_list.get(),
4182 : : uv[ii*2], uv[ii*2+1] ) != CUBIT_SUCCESS)
4183 : : {
4184 : : td_fbp->add_surf_facets( facet_list );
4185 : : td_fbp->set_uv( facet_list.get(), uv[ii*2], uv[ii*2+1] );
4186 : : }
4187 : : }
4188 : : }
4189 : : else
4190 : : {
4191 : : pt->set_uv( uv[ii*2], uv[ii*2+1] );
4192 : : }
4193 : : }
4194 : : isParameterized = CUBIT_TRUE;
4195 : : PRINT_INFO("Surface Parameterization succeeded\n");
4196 : : delete [] sizes;
4197 : : delete [] uv;
4198 : : }
4199 : :
4200 : : // clean up
4201 : :
4202 [ # # ]: 0 : delete [] points;
4203 [ # # ]: 0 : delete [] facets;
4204 : 0 : return CUBIT_SUCCESS;
4205 : : }
4206 : :
4207 : : //===========================================================================
4208 : : //Function Name: compute_curve_tangent
4209 : : //
4210 : : //Member Type: PRIVATE
4211 : : //Descriptoin: compute the tangents to the endpoints of a boundary edge.
4212 : : //===========================================================================
4213 : 0 : CubitStatus FacetEvalTool::compute_curve_tangent(
4214 : : CubitFacetEdge *edge,
4215 : : double min_dot,
4216 : : CubitVector &T0,
4217 : : CubitVector &T3 )
4218 : : {
4219 : :
4220 : 0 : CubitPoint *p0 = edge->point( 0 );
4221 : 0 : CubitPoint *p1 = edge->point( 1 );
4222 : 0 : CubitFacetEdge *prev_edge = next_boundary_edge( edge, p0 );
4223 [ # # ]: 0 : if (prev_edge == NULL) // could be end of a hard line
4224 : : {
4225 [ # # ][ # # ]: 0 : T0 = p1->coordinates() - p0->coordinates();
[ # # ]
4226 : 0 : T0.normalize();
4227 : : }
4228 : : else
4229 : : {
4230 [ # # ]: 0 : CubitPoint *p2 = prev_edge->other_point( p0 );
4231 [ # # ][ # # ]: 0 : CubitVector e0 = p0->coordinates() - p2->coordinates();
[ # # ]
4232 [ # # ][ # # ]: 0 : CubitVector e1 = p1->coordinates() - p0->coordinates();
[ # # ]
4233 [ # # ]: 0 : e0.normalize();
4234 [ # # ]: 0 : e1.normalize();
4235 [ # # ][ # # ]: 0 : if (e0 % e1 >= min_dot)
4236 : : {
4237 [ # # ][ # # ]: 0 : T0 = (p0->coordinates() - p2->coordinates()) +
[ # # ][ # # ]
4238 [ # # ][ # # ]: 0 : (p1->coordinates() - p0->coordinates());
[ # # ][ # # ]
4239 [ # # ]: 0 : T0.normalize();
4240 : : }
4241 : : else
4242 : : {
4243 [ # # ]: 0 : T0 = e1;
4244 : : }
4245 : : }
4246 : :
4247 : :
4248 : 0 : CubitFacetEdge *next_edge = next_boundary_edge( edge, p1 );
4249 [ # # ]: 0 : if (next_edge == NULL) // could be end of a hard line
4250 : : {
4251 [ # # ][ # # ]: 0 : T3 = p1->coordinates() - p0->coordinates();
[ # # ]
4252 : 0 : T3.normalize();
4253 : : }
4254 : : else
4255 : : {
4256 [ # # ]: 0 : CubitPoint *p2 = next_edge->other_point( p1 );
4257 [ # # ][ # # ]: 0 : CubitVector e0 = p2->coordinates() - p1->coordinates();
[ # # ]
4258 [ # # ][ # # ]: 0 : CubitVector e1 = p1->coordinates() - p0->coordinates();
[ # # ]
4259 [ # # ]: 0 : e0.normalize();
4260 [ # # ]: 0 : e1.normalize();
4261 [ # # ][ # # ]: 0 : if (e0 % e1 >= min_dot)
4262 : : {
4263 [ # # ][ # # ]: 0 : T3 = (p2->coordinates() - p1->coordinates()) +
[ # # ][ # # ]
4264 [ # # ][ # # ]: 0 : (p1->coordinates() - p0->coordinates());
[ # # ][ # # ]
4265 [ # # ]: 0 : T3.normalize();
4266 : : }
4267 : : else
4268 : : {
4269 [ # # ]: 0 : T3 = e1;
4270 : : }
4271 : : }
4272 : :
4273 : 0 : return CUBIT_SUCCESS;
4274 : : }
4275 : :
4276 : :
4277 : : //===========================================================================
4278 : : //Function Name: next_boundary_edge
4279 : : //
4280 : : //Member Type: PRIVATE
4281 : : //Descriptoin: given a facet boundary edge and one of its nodes, find the
4282 : : // next edge on the same surface
4283 : : //===========================================================================
4284 : 0 : CubitFacetEdge *FacetEvalTool::next_boundary_edge(
4285 : : CubitFacetEdge *this_edge,
4286 : : CubitPoint *p0 )
4287 : : {
4288 : 0 : CubitFacetEdge *next_edge = NULL;
4289 : :
4290 [ # # ]: 0 : DLIList<CubitFacetEdge*> edge_list;
4291 [ # # ]: 0 : p0->edges( edge_list );
4292 : : int ii;
4293 : :
4294 : 0 : CubitFacetEdge *edge_ptr = NULL;
4295 [ # # ][ # # ]: 0 : for (ii=0; ii<edge_list.size() && next_edge == NULL; ii++)
[ # # ][ # # ]
4296 : : {
4297 [ # # ]: 0 : edge_ptr = edge_list.get_and_step();
4298 [ # # ]: 0 : if (edge_ptr != this_edge)
4299 : : {
4300 [ # # ][ # # ]: 0 : if (edge_ptr->num_adj_facets() <= 1)
4301 : : {
4302 : 0 : next_edge = edge_ptr;
4303 : : }
4304 : : }
4305 : : }
4306 : :
4307 [ # # ]: 0 : return next_edge;
4308 : : }
4309 : :
4310 : : //===========================================================================
4311 : : //Function Name: save
4312 : : //
4313 : : //Member Type: PRIVATE
4314 : : //Description: save the facet eval tool to a cubit file
4315 : : // Assumption: contained facets have been previuosly saved. This function
4316 : : // saves only the facet ids.
4317 : : //===========================================================================
4318 : 275 : CubitStatus FacetEvalTool::save(
4319 : : FILE *fp )
4320 : : {
4321 [ + - ]: 275 : NCubitFile::CIOWrapper cio(fp);
4322 : : typedef NCubitFile::UnsignedInt32 int32;
4323 : :
4324 [ + - ]: 275 : cio.Write(reinterpret_cast<int32*>(&interpOrder), 1);
4325 [ + - ]: 275 : cio.Write(reinterpret_cast<int32*>(&isFlat), 1);
4326 [ - + ]: 275 : int32 is_param = isParameterized ? 1 : 0;
4327 [ + - ]: 275 : cio.Write(&is_param, 1);
4328 : :
4329 [ + - ]: 275 : cio.Write(&myArea, 1);
4330 [ + - ]: 275 : cio.Write(&minDot, 1);
4331 : :
4332 : : // write ids of facets that belong to this tool
4333 : : int32 ii;
4334 : : CubitFacet *facet_ptr;
4335 [ + - ]: 275 : int32 nfacets = myFacetList.size();
4336 [ + - ][ + - ]: 275 : int32* facet_id = new int32 [nfacets];
4337 [ + - ]: 275 : myFacetList.reset();
4338 [ + + ]: 1067 : for (ii=0; ii<nfacets; ii++)
4339 : : {
4340 [ + - ]: 792 : facet_ptr = myFacetList.get_and_step();
4341 [ + - ]: 792 : facet_id[ii] = facet_ptr->id();
4342 : : }
4343 [ + - ]: 275 : cio.Write(&nfacets, 1);
4344 [ + - ]: 275 : if (nfacets > 0)
4345 : : {
4346 [ + - ]: 275 : cio.Write(facet_id, nfacets);
4347 : : }
4348 [ + - ]: 275 : delete [] facet_id;
4349 : :
4350 : : // write ids of edges that belong to this tool
4351 : : CubitFacetEdge *edge_ptr;
4352 [ + - ]: 275 : int32 nedges = myEdgeList.size();
4353 [ + - ][ + - ]: 275 : int32* edge_id = new int32 [nedges];
4354 [ + - ]: 275 : myEdgeList.reset();
4355 [ + + ]: 1969 : for (ii=0; ii<nedges; ii++)
4356 : : {
4357 [ + - ]: 1694 : edge_ptr = myEdgeList.get_and_step();
4358 [ + - ]: 1694 : edge_id[ii] = edge_ptr->id();
4359 : : //volatile int test_int = edge_id[ii] + 1; // this is a test line to look for uninitialized data rjm
4360 : : }
4361 [ + - ]: 275 : cio.Write( &nedges, 1);
4362 [ + - ]: 275 : if (nedges > 0)
4363 : : {
4364 [ + - ]: 275 : cio.Write(edge_id, nedges);
4365 : : }
4366 [ + - ]: 275 : delete [] edge_id;
4367 : :
4368 : : // write ids of points that belong to this tool
4369 : : CubitPoint *point_ptr;
4370 [ + - ]: 275 : int32 npoints = myPointList.size();
4371 [ + - ][ + - ]: 275 : int32* point_id = new int32 [npoints];
4372 [ + - ]: 275 : myPointList.reset();
4373 [ + + ]: 1463 : for (ii=0; ii<npoints; ii++)
4374 : : {
4375 [ + - ]: 1188 : point_ptr = myPointList.get_and_step();
4376 [ + - ]: 1188 : point_id[ii] = point_ptr->id();
4377 : : }
4378 [ + - ]: 275 : cio.Write(&npoints, 1);
4379 [ + - ]: 275 : if (npoints > 0)
4380 : : {
4381 [ + - ]: 275 : cio.Write(point_id, npoints);
4382 : : }
4383 [ + - ]: 275 : delete [] point_id;
4384 : :
4385 [ + - ]: 275 : return CUBIT_SUCCESS;
4386 : : }
4387 : :
4388 : : //===========================================================================
4389 : : //Function Name: restore
4390 : : //
4391 : : //Member Type: PRIVATE
4392 : : //Description: restore a facetevaltool from a CUB file
4393 : : //===========================================================================
4394 : 132 : CubitStatus FacetEvalTool::restore(
4395 : : FILE *fp,
4396 : : unsigned int endian,
4397 : : int num_facets,
4398 : : int num_edges,
4399 : : int num_points,
4400 : : CubitFacet **facets,
4401 : : CubitFacetEdge **edges,
4402 : : CubitPoint **points )
4403 : : {
4404 [ + - ]: 132 : NCubitFile::CIOWrapper cio(endian, fp);
4405 : : typedef NCubitFile::UnsignedInt32 int32;
4406 : :
4407 : : // read interpOrder, isFlat, isParameterized
4408 : : int int_data[3];
4409 [ + - ]: 132 : cio.Read(reinterpret_cast<int32*>(int_data), 3);
4410 : 132 : interpOrder = int_data[0];
4411 : 132 : isFlat = int_data[1];
4412 : 132 : isParameterized = (int_data[2] != 0);
4413 : :
4414 : : // read myArea, minDot
4415 : : double double_data[2];
4416 [ + - ]: 132 : cio.Read( double_data, 2);
4417 : 132 : myArea = double_data[0];
4418 : 132 : minDot = double_data[1];
4419 : :
4420 : 132 : lastFacet = NULL;
4421 : :
4422 : : // read the facet ids and assign facets to this eval tool
4423 : : int nfacets;
4424 [ + - ]: 132 : cio.Read(reinterpret_cast<int32*>(&nfacets), 1);
4425 : 132 : int32 *facet_id = NULL;
4426 [ + - ]: 132 : if (nfacets > 0)
4427 : : {
4428 [ + - ][ + - ]: 132 : facet_id = new int32 [nfacets];
4429 [ + - ]: 132 : cio.Read(facet_id, nfacets);
4430 : : int ii, id;
4431 [ + + ]: 660 : for(ii=0; ii<nfacets; ii++)
4432 : : {
4433 : 528 : id = facet_id[ii];
4434 [ + - ][ - + ]: 528 : if (ii < 0 || ii >= num_facets)
4435 : : {
4436 [ # # ]: 0 : delete [] facet_id;
4437 : 0 : return CUBIT_FAILURE;
4438 : : }
4439 [ + - ]: 528 : myFacetList.append( facets[id] );
4440 [ + - ]: 528 : facets[id]->set_tool_id( toolID );
4441 : : }
4442 [ + - ]: 132 : delete [] facet_id;
4443 : : }
4444 : :
4445 : : // read the edges
4446 : : int nedges;
4447 [ + - ]: 132 : cio.Read(reinterpret_cast<int32*>(&nedges), 1);
4448 : 132 : int32 *edge_id = NULL;
4449 [ + - ]: 132 : if (nedges > 0)
4450 : : {
4451 [ + - ][ + - ]: 132 : edge_id = new int32 [nedges];
4452 [ + - ]: 132 : cio.Read(edge_id, nedges);
4453 : : int ii, id;
4454 [ + + ]: 1188 : for(ii=0; ii<nedges; ii++)
4455 : : {
4456 : 1056 : id = edge_id[ii];
4457 [ + - ][ - + ]: 1056 : if (ii < 0 || ii >= num_edges)
4458 : : {
4459 [ # # ]: 0 : delete [] edge_id;
4460 : 0 : return CUBIT_FAILURE;
4461 : : }
4462 [ + - ]: 1056 : myEdgeList.append( edges[id] );
4463 : : }
4464 [ + - ]: 132 : delete [] edge_id;
4465 : : }
4466 : :
4467 : : // read the points
4468 : : int npoints;
4469 [ + - ]: 132 : cio.Read(reinterpret_cast<int32*>(&npoints), 1);
4470 : 132 : int32 *point_id = NULL;
4471 [ + - ]: 132 : if (npoints > 0)
4472 : : {
4473 [ + - ][ + - ]: 132 : point_id = new int32 [npoints];
4474 [ + - ]: 132 : cio.Read(point_id, npoints);
4475 : : int ii, id;
4476 [ + + ]: 792 : for(ii=0; ii<npoints; ii++)
4477 : : {
4478 : 660 : id = point_id[ii];
4479 [ + - ][ - + ]: 660 : if (ii < 0 || ii >= num_points)
4480 : : {
4481 [ # # ]: 0 : delete [] point_id;
4482 : 0 : return CUBIT_FAILURE;
4483 : : }
4484 [ + - ]: 660 : myPointList.append( points[id] );
4485 : : }
4486 [ + - ]: 132 : delete [] point_id;
4487 : : }
4488 : :
4489 [ + - ][ + - ]: 132 : bounding_box();
4490 [ + - ][ + - ]: 132 : double diag = sqrt(FacetEvalToolUtils::sqr(myBBox->x_range()) +
4491 [ + - ][ + - ]: 132 : FacetEvalToolUtils::sqr(myBBox->y_range()) +
4492 [ + - ][ + - ]: 132 : FacetEvalToolUtils::sqr(myBBox->z_range()));
4493 : 132 : compareTol = 1.0e-3 * diag;
4494 : :
4495 [ + - ]: 132 : return CUBIT_SUCCESS;
4496 : : }
4497 : :
4498 : 0 : CubitStatus FacetEvalTool::get_intersections(CubitVector point1,
4499 : : CubitVector point2,
4500 : : DLIList<CubitVector*>& intersection_list,
4501 : : bool bounded)
4502 : : {
4503 : : //Find the points were the line intersects the bounding box.
4504 [ # # ]: 0 : CubitVector intersect_1;
4505 [ # # ]: 0 : CubitVector intersect_2;
4506 [ # # ]: 0 : CubitBox bbox = *myBBox;
4507 : :
4508 : : //Increase the size of the box in each direction.
4509 : : //Don't use scale because the box may be too thin (planar surface).
4510 [ # # ][ # # ]: 0 : double offset = 2.0 * (point1 - point2).length();
4511 : :
4512 [ # # ]: 0 : CubitVector min;
4513 [ # # ][ # # ]: 0 : min.x( bbox.min_x() - offset );
4514 [ # # ][ # # ]: 0 : min.y( bbox.min_y() - offset );
4515 [ # # ][ # # ]: 0 : min.z( bbox.min_z() - offset );
4516 [ # # ]: 0 : CubitVector max;
4517 [ # # ][ # # ]: 0 : max.x( bbox.max_x() + offset );
4518 [ # # ][ # # ]: 0 : max.y( bbox.max_y() + offset );
4519 [ # # ][ # # ]: 0 : max.z( bbox.max_z() + offset );
4520 : :
4521 [ # # ]: 0 : bbox.reset( min, max );
4522 : : int box_intersections = FacetDataUtil::get_bbox_intersections( point1, point2, bbox,
4523 [ # # ]: 0 : intersect_1, intersect_2 );
4524 : :
4525 : : //The bounding box is larger than the surface we are checking.
4526 : : //This means that if there are less than two intersections
4527 : : //the line will not intersect the surface.
4528 [ # # ]: 0 : if( 2 > box_intersections )
4529 : 0 : return CUBIT_SUCCESS;
4530 : :
4531 [ # # ]: 0 : bbox.reset( intersect_1 );
4532 [ # # ]: 0 : bbox |= intersect_2;
4533 : :
4534 : : //Find the facets that are intersected by the bbox that was just created.
4535 [ # # ][ # # ]: 0 : DLIList<CubitFacet*> search_facets;
4536 [ # # ]: 0 : if( aTree )
4537 : : {
4538 : : //Get the facets from the tree.
4539 [ # # ]: 0 : aTree->find( bbox, search_facets );
4540 : : }
4541 : : else
4542 [ # # ]: 0 : search_facets = myFacetList;
4543 : :
4544 : : int ii;
4545 [ # # ][ # # ]: 0 : for( ii = search_facets.size(); ii > 0; ii-- )
4546 : : {
4547 [ # # ]: 0 : CubitFacet* test_facet = search_facets.get_and_step();
4548 : :
4549 [ # # ]: 0 : CubitVector intersection;
4550 [ # # ]: 0 : CubitVector area_coord;
4551 : : CubitPointContainment contain = FacetDataUtil::intersect_facet( intersect_1, intersect_2, test_facet,
4552 [ # # ]: 0 : intersection, area_coord );
4553 : :
4554 [ # # ]: 0 : if( bounded )
4555 : : {
4556 [ # # ]: 0 : CubitVector dir1 = point2 - point1;
4557 [ # # ]: 0 : CubitVector dir2 = intersection - point1;
4558 : :
4559 [ # # ][ # # ]: 0 : if( dir2.length_squared() > (GEOMETRY_RESABS * GEOMETRY_RESABS) )
4560 : : {
4561 [ # # ][ # # ]: 0 : if( dir1 % dir2 < 0 ||
[ # # ][ # # ]
4562 [ # # ][ # # ]: 0 : ( ( dir2.length_squared() - dir1.length_squared() ) >
4563 : : ( GEOMETRY_RESABS * GEOMETRY_RESABS ) ) )
4564 : : {
4565 : : //The inserction point is not between the two end points.
4566 : 0 : contain = CUBIT_PNT_OUTSIDE;
4567 : : }
4568 : : }
4569 : : }
4570 : :
4571 [ # # ][ # # ]: 0 : if( CUBIT_PNT_BOUNDARY == contain ||
4572 : : CUBIT_PNT_INSIDE == contain )
4573 : : {
4574 : : //The point intersects the facets.
4575 [ # # ][ # # ]: 0 : CubitVector* new_intersection = new CubitVector;
4576 [ # # ]: 0 : *new_intersection = intersection;
4577 [ # # ]: 0 : intersection_list.append( new_intersection );
4578 : : }
4579 : : }
4580 : :
4581 : : //Remove duplicate intersections.
4582 [ # # ]: 0 : intersection_list.reset();
4583 [ # # ][ # # ]: 0 : for( ii = 0; ii < intersection_list.size(); ii++ )
4584 : : {
4585 [ # # ]: 0 : CubitVector* base_vec = intersection_list.next(ii);
4586 [ # # ]: 0 : if( NULL == base_vec )
4587 : 0 : continue;
4588 : :
4589 : : int jj;
4590 [ # # ][ # # ]: 0 : for( jj = ii+1; jj < intersection_list.size(); jj++ )
4591 : : {
4592 [ # # ]: 0 : CubitVector* compare_vec = intersection_list.next(jj);
4593 [ # # ]: 0 : if( NULL != compare_vec )
4594 : : {
4595 [ # # ][ # # ]: 0 : if( base_vec->distance_between_squared( *compare_vec ) < GEOMETRY_RESABS * GEOMETRY_RESABS )
4596 : : {
4597 [ # # ]: 0 : intersection_list.step(jj);
4598 [ # # ]: 0 : delete intersection_list.get();
4599 [ # # ]: 0 : intersection_list.change_to( NULL );
4600 [ # # ]: 0 : intersection_list.reset();
4601 : : }
4602 : : }
4603 : : }
4604 : : }
4605 [ # # ]: 0 : intersection_list.remove_all_with_value( NULL );
4606 : :
4607 [ # # ]: 0 : return CUBIT_SUCCESS;
4608 : : }
4609 : :
4610 : 0 : int FacetEvalTool::intersect_ray( CubitVector &origin, CubitVector &direction, CubitFacet* facet, CubitVector* point, double &distance )
4611 : : {
4612 : : // This algorithm can be found at http://geometryalgorithms.com/
4613 : :
4614 [ # # ]: 0 : CubitVector n; // triangle vectors
4615 [ # # ][ # # ]: 0 : CubitVector w0, w; // ray vectors
4616 : : double a, b; // params to calc ray-plane intersect
4617 : :
4618 : : // get triangle edge vectors and plane normal
4619 : : //CubitVector normal = facet->normal();
4620 [ # # ]: 0 : CubitPoint *pt0 = facet->point(0);
4621 [ # # ]: 0 : CubitPoint *pt1 = facet->point(1);
4622 [ # # ]: 0 : CubitPoint *pt2 = facet->point(2);
4623 : 0 : double tol = GEOMETRY_RESABS;
4624 : :
4625 [ # # ][ # # ]: 0 : CubitVector u( pt1->x() - pt0->x(),
4626 [ # # ][ # # ]: 0 : pt1->y() - pt0->y(),
4627 [ # # ][ # # ]: 0 : pt1->z() - pt0->z()); //(*p1-*p0);
[ # # ]
4628 [ # # ][ # # ]: 0 : CubitVector v( pt2->x() - pt0->x(),
4629 [ # # ][ # # ]: 0 : pt2->y() - pt0->y(),
4630 [ # # ][ # # ]: 0 : pt2->z() - pt0->z()); // = (*p2-*p0);
[ # # ]
4631 : :
4632 : : //u = T.V1 - T.V0;
4633 : : //v = T.V2 - T.V0;
4634 [ # # ][ # # ]: 0 : n = u * v; // cross product
4635 [ # # ][ # # ]: 0 : if (n.length_squared() == 0) // triangle is degenerate
4636 : 0 : return -1; // do not deal with this case
4637 : :
4638 : : //dir = R.P1 - R.P0; // ray direction vector
4639 : : //w0 = R.P0 - T.V0;
4640 [ # # ][ # # ]: 0 : w0 = CubitVector(origin.x() - pt0->x(),
[ # # ]
4641 [ # # ][ # # ]: 0 : origin.y() - pt0->y(),
4642 [ # # ][ # # ]: 0 : origin.z() - pt0->z());
[ # # ]
4643 : :
4644 [ # # ]: 0 : a = -(n%w0);
4645 [ # # ]: 0 : b = (n%direction);
4646 [ # # ]: 0 : if (fabs(b) < tol) { // ray is parallel to triangle plane
4647 [ # # ]: 0 : if (a == 0) // ray lies in triangle plane
4648 : 0 : return 2;
4649 : 0 : else return 0; // ray disjoint from plane
4650 : : }
4651 : :
4652 : : // get intersect point of ray with triangle plane
4653 : 0 : distance = a / b;
4654 [ # # ]: 0 : if (distance < 0.0) // ray goes away from triangle
4655 : 0 : return 0; // => no intersect
4656 : : // for a segment, also test if (r > 1.0) => no intersect
4657 : :
4658 [ # # ][ # # ]: 0 : point->set(origin + distance * direction); // intersect point of ray and plane
[ # # ]
4659 : :
4660 : : // the distance we want to return is real distance, not distance/magnitude
4661 [ # # ]: 0 : distance *= direction.length();
4662 : :
4663 : : // is point inside facet?
4664 : : double uu, uv, vv, wu, wv, D;
4665 [ # # ]: 0 : uu = u%u;
4666 [ # # ]: 0 : uv = u%v;
4667 [ # # ]: 0 : vv = v%v;
4668 : : //w = *I - T.V0;
4669 [ # # ][ # # ]: 0 : w = CubitVector(point->x() - pt0->x(),
[ # # ]
4670 [ # # ][ # # ]: 0 : point->y() - pt0->y(),
4671 [ # # ][ # # ]: 0 : point->z() - pt0->z());
[ # # ]
4672 [ # # ]: 0 : wu = w%u;
4673 [ # # ]: 0 : wv = w%v;
4674 : 0 : D = uv * uv - uu * vv;
4675 : :
4676 : : // get and test parametric coords
4677 : : double s, t;
4678 : 0 : s = (uv * wv - vv * wu) / D;
4679 [ # # ][ # # ]: 0 : if (s < 0.0 || s > 1.0) // point is outside facet
4680 : 0 : return 0;
4681 : 0 : t = (uv * wu - uu * wv) / D;
4682 [ # # ][ # # ]: 0 : if (t < 0.0 || (s + t) > 1.0) // point is outside facet
4683 : 0 : return 0;
4684 : :
4685 : 0 : return 1; // point is in facet
4686 : :
4687 : : }
4688 : :
4689 : 0 : int FacetEvalTool::intersect_ray( CubitVector &origin, CubitVector &direction, CubitFacetEdge* facet_edge, CubitVector* point, double &hit_distance )
4690 : : {
4691 : : // This algorithm can be found at http://geometryalgorithms.com/
4692 : 0 : double tol = GEOMETRY_RESABS;
4693 : :
4694 [ # # ]: 0 : CubitPoint* p0 = facet_edge->point(0);
4695 [ # # ]: 0 : CubitPoint* p1 = facet_edge->point(1);
4696 [ # # ][ # # ]: 0 : CubitVector u0 = CubitVector(p0->x(), p0->y(), p0->z());
[ # # ][ # # ]
4697 [ # # ][ # # ]: 0 : CubitVector u1 = CubitVector(p1->x(), p1->y(), p1->z());
[ # # ][ # # ]
4698 : :
4699 [ # # ]: 0 : CubitVector u = CubitVector(u0, u1);
4700 [ # # ]: 0 : CubitVector v = direction;
4701 [ # # ]: 0 : v.normalize();
4702 : :
4703 [ # # ]: 0 : CubitVector w = CubitVector(origin, u0);
4704 : :
4705 : : double sc, tc; // sc is fraction along facet edge, tc is distance along ray
4706 : :
4707 [ # # ]: 0 : double a = u%u; // always >= 0
4708 [ # # ]: 0 : double b = u%v;
4709 [ # # ]: 0 : double c = v%v; // always >= 0
4710 [ # # ]: 0 : double d = u%w;
4711 [ # # ]: 0 : double e = v%w;
4712 : 0 : double D = a*c - b*b; // always >= 0
4713 : :
4714 : : // compute the line parameters of the two closest points
4715 [ # # ]: 0 : if (D < tol)
4716 : : {
4717 : : // the lines are almost parallel
4718 : 0 : sc = 0.0;
4719 [ # # ]: 0 : tc = (b>c ? d/b : e/c); // use the largest denominator
4720 : : }
4721 : : else
4722 : : {
4723 : 0 : sc = (b*e - c*d) / D;
4724 : 0 : tc = (a*e - b*d) / D;
4725 : : }
4726 : :
4727 : : // get the difference of the two closest points
4728 [ # # ][ # # ]: 0 : CubitVector dP = CubitVector(w + (sc * u) - (tc * v)); // = <0 0 0> if intersection
[ # # ][ # # ]
4729 : :
4730 [ # # ]: 0 : double distance = sqrt(dP % dP); // return the closest distance (0 if intersection)
4731 : :
4732 [ # # ][ # # ]: 0 : point->set(u0 + (sc * u));
[ # # ]
4733 : 0 : hit_distance = tc; //distance from origin to intersection point
4734 : :
4735 [ # # ]: 0 : if (distance < tol)
4736 : : {
4737 : : //check if parallel (infinite intersection)
4738 [ # # ]: 0 : if (D < tol)
4739 : 0 : return 2;
4740 : : //check if on edge
4741 [ # # ][ # # ]: 0 : if (sc <= 1.0 && sc >= 0.0)
4742 : 0 : return 1;
4743 : : else
4744 : 0 : return 0;
4745 : : }
4746 : :
4747 : 0 : return 0;
4748 : : }
4749 : :
4750 : 11 : double FacetEvalTool::contained_volume( DLIList<CubitFacet *> &facets )
4751 : : {
4752 [ + - ][ + - ]: 11 : CubitVector p0, p1, p2, p3, normal;
[ + - ][ + - ]
[ + - ]
4753 : 11 : double volume = 0.0;
4754 : :
4755 [ + - ]: 11 : facets.reset();
4756 [ + - ][ + + ]: 143 : for (int i = facets.size(); i--; )
4757 : : {
4758 [ + - ]: 132 : CubitFacet* facet = facets.get_and_step();
4759 [ + - ][ + - ]: 132 : p1 = facet->point(0)->coordinates();
[ + - ]
4760 [ + - ][ + - ]: 132 : p2 = facet->point(1)->coordinates();
[ + - ]
4761 [ + - ][ + - ]: 132 : p3 = facet->point(2)->coordinates();
[ + - ]
4762 [ + - ][ + - ]: 132 : normal = (p3 - p1) * (p2 - p1);
[ + - ][ + - ]
4763 : :
4764 [ + - ]: 132 : double two_area = normal.length();
4765 [ + - ]: 132 : if (two_area > CUBIT_RESABS )
4766 : : {
4767 [ + - ]: 132 : normal /= two_area;
4768 : :
4769 [ + - ][ + - ]: 132 : double height = normal % (p0 - p1);
4770 : 132 : double vol = two_area * height;
4771 : 132 : volume += vol;
4772 : : }
4773 : : }
4774 : :
4775 : 11 : volume /= 6.0;
4776 : 11 : return volume;
4777 [ + - ][ + - ]: 6540 : }
4778 : :
|
