Branch data Line data Source code
1 : : //=======================================================================
2 : : //
3 : : // File: CubitFacetEdge
4 : : // Description: used for edges of CubitFacets. These are optional
5 : : // and used only if specific information must be stored at
6 : : // the edge common to two facets.
7 : : // Notes: See notes in CubitFacetEdge.hpp
8 : : // Owner: sjowen
9 : : //
10 : : //=======================================================================
11 : :
12 : : #include "CubitFacetEdge.hpp"
13 : : #include "CubitFacet.hpp"
14 : : #include "CubitPoint.hpp"
15 : : #include "CubitVector.hpp"
16 : : #include "GeometryDefines.h"
17 : : #include "FacetEvalTool.hpp"
18 : : #include "GfxDebug.hpp"
19 : : #include "IntersectionTool.hpp"
20 : : #include "CubitMessage.hpp"
21 : :
22 : :
23 : : //======================================================================
24 : : // Function: CubitFacetEdge (PUBLIC)
25 : : // Description: constructor
26 : : // Author: sjowen
27 : : // Date: 8/00
28 : : //======================================================================
29 : 2486 : CubitFacetEdge::CubitFacetEdge( )
30 [ + - ][ + + ]: 9944 : : bezierOrder(0), markedFlag(0), isFeature(0), isFlipped(CUBIT_FALSE)
31 : : {
32 : 2486 : }
33 : :
34 : :
35 : : //======================================================================
36 : : // Function: ~CubitFacetEdge (PUBLIC)
37 : : // Description: destructor
38 : : // Author: sjowen
39 : : // Date: 8/00
40 : : //======================================================================
41 : 0 : CubitFacetEdge::~CubitFacetEdge()
42 : : {
43 [ # # ]: 0 : }
44 : :
45 : : //======================================================================
46 : : // Function: set_control_points (PUBLIC)
47 : : // Description: set the control points from an array of doubles
48 : : // Author: sjowen
49 : : // Date: 8/00
50 : : //======================================================================
51 : 792 : void CubitFacetEdge::set_control_points( const double *ctrl_pt_array )
52 : : {
53 : : int ii;
54 [ + + ]: 3168 : for (ii=0; ii<3; ii++)
55 : : {
56 : 2376 : controlPoints[ii].x( ctrl_pt_array[ii*3] );
57 : 2376 : controlPoints[ii].y( ctrl_pt_array[ii*3+1] );
58 : 2376 : controlPoints[ii].z( ctrl_pt_array[ii*3+2] );
59 : : }
60 : 792 : }
61 : :
62 : : //======================================================================
63 : : // Function: control_points (PUBLIC)
64 : : // Description: return the Bezier control points (including the end points)
65 : : // The order of the bezier is returned.
66 : : // Author: sjowen
67 : : // Date: 8/00
68 : : //======================================================================
69 : 330 : int CubitFacetEdge::control_points( CubitVector *ctrl_pts )
70 : : {
71 [ + - ]: 330 : DLIList<CubitPoint*> my_points;
72 [ + - ]: 330 : points(my_points);
73 : :
74 [ + - ][ + - ]: 330 : ctrl_pts[0] = my_points.get()->coordinates();
[ + - ]
75 [ + + ]: 1320 : for(int i=0; i<bezierOrder-1; i++) {
76 [ + - ]: 990 : ctrl_pts[i+1] = controlPoints[i];
77 : : }
78 [ + - ][ + - ]: 330 : ctrl_pts[bezierOrder] = my_points.next()->coordinates();
[ + - ]
79 [ + - ]: 330 : return bezierOrder;
80 : : }
81 : :
82 : : //======================================================================
83 : : // Function: control_points (PUBLIC)
84 : : // Description: return the control points on the edge of a facet based
85 : : // on its edge use direction
86 : : // Author: sjowen
87 : : // Date: 8/00
88 : : //======================================================================
89 : 1584 : CubitStatus CubitFacetEdge::control_points(
90 : : CubitFacet *facet, CubitVector *ctrl_pts )
91 : : {
92 : 1584 : int index = -1;
93 : 1584 : CubitBoolean found = CUBIT_FALSE;
94 [ + + ][ + + ]: 4752 : for (int i=0; i<3 && !found; i++) {
95 [ + - ][ + + ]: 3168 : if (this == facet->edge(i)) {
96 : 1584 : index = i;
97 : 1584 : found = CUBIT_TRUE;
98 : : }
99 : : }
100 [ - + ]: 1584 : if (!found) {
101 : 0 : return CUBIT_FAILURE;
102 : : }
103 : :
104 [ + - ]: 1584 : DLIList<CubitPoint*> my_points;
105 [ + - ]: 1584 : points(my_points);
106 : :
107 [ + - ]: 1584 : switch (facet->edge_use(index)) {
[ + + - ]
108 : : case 1:
109 [ + - ][ + - ]: 792 : ctrl_pts[0] = my_points.get()->coordinates();
[ + - ]
110 [ + - ]: 792 : ctrl_pts[1] = controlPoints[0];
111 [ + - ]: 792 : ctrl_pts[2] = controlPoints[1];
112 [ + - ]: 792 : ctrl_pts[3] = controlPoints[2];
113 [ + - ][ + - ]: 792 : ctrl_pts[4] = my_points.next()->coordinates();
[ + - ]
114 : 792 : break;
115 : : case -1:
116 [ + - ][ + - ]: 792 : ctrl_pts[0] = my_points.next()->coordinates();
[ + - ]
117 [ + - ]: 792 : ctrl_pts[1] = controlPoints[2];
118 [ + - ]: 792 : ctrl_pts[2] = controlPoints[1];
119 [ + - ]: 792 : ctrl_pts[3] = controlPoints[0];
120 [ + - ][ + - ]: 792 : ctrl_pts[4] = my_points.get()->coordinates();
[ + - ]
121 : 792 : break;
122 : : default:
123 : 0 : return CUBIT_FAILURE;
124 : : }
125 [ + - ]: 1584 : return CUBIT_SUCCESS;
126 : : }
127 : :
128 : : //===========================================================================
129 : : //Function Name: evaluate_position
130 : : //
131 : : //Member Type: PUBLIC
132 : : //Description: evaluate the facet edge at a position
133 : : // eval_tangent is NULL if tangent not needed
134 : : //===========================================================================
135 : 0 : CubitStatus CubitFacetEdge::evaluate_position( const CubitVector &start_position,
136 : : CubitVector *eval_point,
137 : : CubitVector *eval_tangent)
138 : : {
139 : 0 : CubitStatus stat = CUBIT_SUCCESS;
140 : :
141 : : // find the adjacent facet
142 : :
143 : 0 : CubitFacet *facet_ptr = this->adj_facet( 0 );
144 : :
145 : : // If there is none or this is a linear representation -
146 : : // then project to the linear edge
147 : :
148 [ # # ][ # # ]: 0 : if (!facet_ptr || facet_ptr->eval_order() == 0 || facet_ptr->is_flat())
[ # # ][ # # ]
149 : : {
150 [ # # ]: 0 : if (eval_point)
151 : : {
152 : 0 : closest_point(start_position, *eval_point);
153 : : }
154 [ # # ]: 0 : if (eval_tangent)
155 : : {
156 [ # # ][ # # ]: 0 : *eval_tangent = point(1)->coordinates() - point(0)->coordinates();
[ # # ][ # # ]
157 : 0 : (*eval_tangent).normalize();
158 : : }
159 : : }
160 : : else
161 : : {
162 [ # # ][ # # ]: 0 : int vert0 = facet_ptr->point_index( point(0) );
163 [ # # ][ # # ]: 0 : int vert1 = facet_ptr->point_index( point(1) );
164 [ # # ][ # # ]: 0 : CubitVector pt_on_plane, close_point;
165 [ # # ]: 0 : CubitVector start = start_position;
166 : : double dist_to_plane;
167 : : CubitBoolean outside_facet;
168 : : FacetEvalTool::project_to_facet_plane( facet_ptr, start,
169 [ # # ]: 0 : pt_on_plane, dist_to_plane );
170 : : stat = FacetEvalTool::project_to_facetedge( facet_ptr,
171 : : vert0, vert1,
172 : : start,
173 : : pt_on_plane,
174 : : close_point,
175 [ # # ]: 0 : outside_facet );
176 [ # # ]: 0 : if (eval_point)
177 : : {
178 [ # # ]: 0 : *eval_point = close_point;
179 : : }
180 [ # # ]: 0 : if (eval_tangent)
181 : : {
182 [ # # ][ # # ]: 0 : CubitVector edvec = point(1)->coordinates() - point(0)->coordinates();
[ # # ][ # # ]
[ # # ]
183 [ # # ]: 0 : edvec.normalize();
184 [ # # ]: 0 : CubitVector areacoord;
185 [ # # ]: 0 : FacetEvalTool::facet_area_coordinate( facet_ptr, close_point, areacoord );
186 [ # # ]: 0 : FacetEvalTool::eval_facet_normal(facet_ptr, areacoord, *eval_tangent);
187 [ # # ]: 0 : CubitVector cross = edvec * *eval_tangent;
188 [ # # ][ # # ]: 0 : *eval_tangent = *eval_tangent * cross;
189 [ # # ]: 0 : (*eval_tangent).normalize();
190 : : }
191 : : }
192 : :
193 : 0 : return stat;
194 : : }
195 : :
196 : : //===========================================================================
197 : : //Function Name: evaluate
198 : : //
199 : : //Member Type: PUBLIC
200 : : //Description: evaluate the facet at area coordinates
201 : : // eval_normal is NULL if normal not needed
202 : : //Note: t is a value from -1 to 1. t=0 is the midpoint
203 : : //===========================================================================
204 : 0 : CubitStatus CubitFacetEdge::evaluate( double &t,
205 : : CubitVector *eval_point,
206 : : CubitVector *eval_tangent )
207 : : {
208 : 0 : CubitStatus stat = CUBIT_SUCCESS;
209 : :
210 : : // project the position to the linear edge
211 : :
212 : 0 : double tt = (t + 1) * 0.5;
213 [ # # ]: 0 : if (tt <= 0.0) tt = 0.0;
214 [ # # ]: 0 : if (tt >= 1.0) tt = 1.0;
215 [ # # ][ # # ]: 0 : *eval_point = point(0)->coordinates() +
[ # # ]
216 [ # # ][ # # ]: 0 : tt * (point(1)->coordinates() - point(0)->coordinates());
[ # # ][ # # ]
[ # # ]
217 : :
218 : : // evaluate the point on the facet (if the order is higher than 0)
219 : :
220 : 0 : CubitFacet *facet_ptr = this->adj_facet( 0 );
221 [ # # ][ # # ]: 0 : if (!facet_ptr || facet_ptr->is_flat())
[ # # ]
222 : : {
223 [ # # ]: 0 : if (eval_tangent)
224 : : {
225 [ # # ][ # # ]: 0 : *eval_tangent = point(1)->coordinates() - point(0)->coordinates();
[ # # ][ # # ]
226 : 0 : (*eval_tangent).normalize();
227 : : }
228 : : }
229 : : else
230 : : {
231 [ # # ]: 0 : CubitVector areacoord;
232 [ # # ]: 0 : FacetEvalTool::facet_area_coordinate( facet_ptr, *eval_point, areacoord );
233 [ # # ]: 0 : stat = facet_ptr->evaluate( areacoord, eval_point, eval_tangent );
234 [ # # ]: 0 : if (stat != CUBIT_SUCCESS)
235 : 0 : return stat;
236 [ # # ]: 0 : if (eval_tangent)
237 : : {
238 [ # # ][ # # ]: 0 : CubitVector edvec = point(1)->coordinates() - point(0)->coordinates();
[ # # ][ # # ]
[ # # ]
239 [ # # ]: 0 : edvec.normalize();
240 [ # # ]: 0 : CubitVector cross = edvec * *eval_tangent;
241 [ # # ][ # # ]: 0 : *eval_tangent = *eval_tangent * cross;
242 [ # # ]: 0 : (*eval_tangent).normalize();
243 : : }
244 : : }
245 : 0 : return stat;
246 : : }
247 : :
248 : : //===========================================================================
249 : : //Function Name: evaluate_single
250 : : //
251 : : //Member Type: PUBLIC
252 : : //Description: evaluate edge not associated with a facet.
253 : : //Note: t is a value from -1 to 1.
254 : : //===========================================================================
255 : 0 : CubitStatus CubitFacetEdge::evaluate_single(double &t,
256 : : CubitVector *outv)
257 : : {
258 [ # # ][ # # ]: 0 : CubitVector P0, P1;
259 : : double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
260 [ # # ]: 0 : DLIList<CubitPoint*> my_points;
261 : : // project the position to the linear edge
262 : :
263 : 0 : double tt = (t + 1) * 0.5;
264 [ # # ]: 0 : if (tt <= 0.0) tt = 0.0;
265 [ # # ]: 0 : if (tt >= 1.0) tt = 1.0;
266 : :
267 [ # # ]: 0 : points(my_points);
268 : :
269 [ # # ][ # # ]: 0 : P0 = my_points.get()->coordinates();
[ # # ]
270 [ # # ][ # # ]: 0 : P1 = my_points.next()->coordinates();
[ # # ]
271 : :
272 : 0 : t2 = tt*tt;
273 : 0 : t3 = t2*tt;
274 : 0 : t4 = t3*tt;
275 : 0 : one_minus_t = 1.-tt;
276 : 0 : one_minus_t2 = one_minus_t*one_minus_t;
277 : 0 : one_minus_t3 = one_minus_t2*one_minus_t;
278 : 0 : one_minus_t4 = one_minus_t3*one_minus_t;
279 : :
280 [ # # ][ # # ]: 0 : *outv = one_minus_t4*P0 +
281 [ # # ][ # # ]: 0 : 4.*one_minus_t3*tt* controlPoints[0] +
282 [ # # ][ # # ]: 0 : 6.*one_minus_t2*t2*controlPoints[1] +
283 [ # # ][ # # ]: 0 : 4.*one_minus_t* t3*controlPoints[2] +
284 [ # # ][ # # ]: 0 : t4*P1;
285 : :
286 [ # # ]: 0 : return CUBIT_SUCCESS;
287 : : }
288 : :
289 : : //===========================================================================
290 : : //Function Name: evaluate_single_tangent
291 : : //
292 : : //Member Type: PUBLIC
293 : : //Description: evaluate tangent to edge not associated with a facet.
294 : : //Note: t is a value from -1 to 1.
295 : : //===========================================================================
296 : 0 : CubitStatus CubitFacetEdge::evaluate_single_tangent(double &t,
297 : : CubitVector *outv)
298 : : {
299 [ # # ][ # # ]: 0 : CubitVector P0, P1;
300 : : double t3, t2, one_minus_t, one_minus_t2, one_minus_t3;
301 [ # # ]: 0 : DLIList<CubitPoint*> my_points;
302 : : // project the position to the linear edge
303 : :
304 : 0 : double tt = (t + 1) * 0.5;
305 [ # # ]: 0 : if (tt <= 0.0) tt = 0.0;
306 [ # # ]: 0 : if (tt >= 1.0) tt = 1.0;
307 : :
308 [ # # ]: 0 : points(my_points);
309 : :
310 [ # # ][ # # ]: 0 : P0 = my_points.get()->coordinates();
[ # # ]
311 [ # # ][ # # ]: 0 : P1 = my_points.next()->coordinates();
[ # # ]
312 : :
313 : 0 : t2 = tt*tt;
314 : 0 : t3 = t2*tt;
315 : 0 : one_minus_t = 1.-tt;
316 : 0 : one_minus_t2 = one_minus_t*one_minus_t;
317 : 0 : one_minus_t3 = one_minus_t2*one_minus_t;
318 : :
319 [ # # ][ # # ]: 0 : *outv = -4.*one_minus_t3*P0 +
320 [ # # ][ # # ]: 0 : 4.*(one_minus_t3 -3.*tt*one_minus_t2)*controlPoints[0] +
321 [ # # ][ # # ]: 0 : 12.*(tt*one_minus_t2 - t2*one_minus_t)*controlPoints[1] +
322 [ # # ][ # # ]: 0 : 4.*(3.*t2*one_minus_t - t3)*controlPoints[2] +
323 [ # # ][ # # ]: 0 : 4.*t3*P1;
324 : :
325 [ # # ]: 0 : return CUBIT_SUCCESS;
326 : : }
327 : :
328 : : //===========================================================================
329 : : //Function Name: evaluate_2nd_derivative
330 : : //
331 : : //Member Type: PUBLIC
332 : : //Description: evaluate tangent to edge not associated with a facet.
333 : : //Note: t is a value from -1 to 1.
334 : : //===========================================================================
335 : 0 : CubitStatus CubitFacetEdge::evaluate_2nd_derivative(double &t,
336 : : CubitVector *outv)
337 : : {
338 [ # # ][ # # ]: 0 : CubitVector P0, P1, second_d;
[ # # ]
339 [ # # ]: 0 : DLIList<CubitPoint*> my_points;
340 : : double val;
341 : : // project the position to the linear edge
342 : :
343 : 0 : double tt = (t + 1) * 0.5;
344 [ # # ]: 0 : if (tt <= 0.0) tt = 0.0;
345 [ # # ]: 0 : if (tt >= 1.0) tt = 1.0;
346 : :
347 [ # # ]: 0 : points(my_points);
348 : :
349 [ # # ][ # # ]: 0 : P0 = my_points.get()->coordinates();
[ # # ]
350 [ # # ][ # # ]: 0 : P1 = my_points.next()->coordinates();
[ # # ]
351 : :
352 [ # # ]: 0 : val = 12.*(1.-tt)*(1.-tt)*P0.x() -
353 [ # # ]: 0 : 24.*(2.*tt*tt - 3.*tt + 1.)*controlPoints[0].x() +
354 [ # # ]: 0 : 12.*(6.*tt*tt - 6.*tt + 1.)*controlPoints[1].x() +
355 [ # # ]: 0 : 24.*(tt - 2.*tt*tt)*controlPoints[2].x() +
356 [ # # ]: 0 : 12.*tt*tt*P1.x();
357 [ # # ]: 0 : second_d.x(val);
358 [ # # ]: 0 : val = 12.*(1.-tt)*(1.-tt)*P0.y() -
359 [ # # ]: 0 : 24.*(2.*tt*tt - 3.*tt + 1.)*controlPoints[0].y() +
360 [ # # ]: 0 : 12.*(6.*tt*tt - 6.*tt + 1.)*controlPoints[1].y() +
361 [ # # ]: 0 : 24.*(tt - 2.*tt*tt)*controlPoints[2].y() +
362 [ # # ]: 0 : 12.*tt*tt*P1.y();
363 [ # # ]: 0 : second_d.y(val);
364 [ # # ]: 0 : val = 12.*(1.-tt)*(1.-tt)*P0.z() -
365 [ # # ]: 0 : 24.*(2.*tt*tt - 3.*tt + 1.)*controlPoints[0].z() +
366 [ # # ]: 0 : 12.*(6.*tt*tt - 6.*tt + 1.)*controlPoints[1].z() +
367 [ # # ]: 0 : 24.*(tt - 2.*tt*tt)*controlPoints[2].z() +
368 [ # # ]: 0 : 12.*tt*tt*P1.z();
369 [ # # ]: 0 : second_d.z(val);
370 [ # # ]: 0 : *outv = second_d;
371 : :
372 [ # # ]: 0 : return CUBIT_SUCCESS;
373 : : }
374 : :
375 : : //===========================================================================
376 : : //Function Name: closest_point
377 : : //
378 : : //Member Type: PUBLIC
379 : : //Description: return the closest point on segment defined by the edge
380 : : //===========================================================================
381 : 0 : CubitStatus CubitFacetEdge::closest_point(const CubitVector &point,
382 : : CubitVector &closest_point )
383 : : {
384 : : //CubitStatus rv = CUBIT_SUCCESS;
385 [ # # ]: 0 : CubitPoint *pt0 = this->point(0);
386 [ # # ]: 0 : CubitPoint *pt1 = this->point(1);
387 : :
388 : : // the edge vector
389 : :
390 [ # # ][ # # ]: 0 : CubitVector e0 ( pt1->x() - pt0->x(),
391 [ # # ][ # # ]: 0 : pt1->y() - pt0->y(),
392 [ # # ][ # # ]: 0 : pt1->z() - pt0->z() );
[ # # ]
393 [ # # ]: 0 : double elen = e0.normalize();
394 : :
395 : : // vector from vert0 to point
396 : :
397 [ # # ][ # # ]: 0 : CubitVector v0 ( point.x() - pt0->x(),
398 [ # # ][ # # ]: 0 : point.y() - pt0->y(),
399 [ # # ][ # # ]: 0 : point.z() - pt0->z() );
[ # # ]
400 : :
401 : : // project to edge
402 : :
403 [ # # ]: 0 : double len = v0%e0;
404 [ # # ]: 0 : if (len <= 0.0)
405 : : {
406 [ # # ][ # # ]: 0 : closest_point = pt0->coordinates();
407 : : }
408 [ # # ]: 0 : else if( len >= elen )
409 : : {
410 [ # # ][ # # ]: 0 : closest_point = pt1->coordinates();
411 : : }
412 : : else
413 : : {
414 [ # # ][ # # ]: 0 : closest_point.x ( pt0->x() + e0.x() * len );
[ # # ]
415 [ # # ][ # # ]: 0 : closest_point.y ( pt0->y() + e0.y() * len );
[ # # ]
416 [ # # ][ # # ]: 0 : closest_point.z ( pt0->z() + e0.z() * len );
[ # # ]
417 : : }
418 : :
419 : 0 : return CUBIT_SUCCESS;
420 : : }
421 : :
422 : : //===========================================================================
423 : : //Function Name: intersect
424 : : //
425 : : //Member Type: PUBLIC
426 : : //Description: intersect the edge with a segment. Assumes segment and edge
427 : : // are on the same plane (project to facet plane first)
428 : : //===========================================================================
429 : 0 : CubitStatus CubitFacetEdge::intersect(
430 : : CubitVector &aa, CubitVector &bb, // end point of segment
431 : : CubitVector &norm, // normal of the common plane
432 : : CubitVector &qq, // return the intersection point
433 : : CubitBoolean &does_intersect ) // return status of intersection
434 : : {
435 : :
436 [ # # ]: 0 : CubitPoint *pt0 = this->point(0);
437 [ # # ]: 0 : CubitPoint *pt1 = this->point(1);
438 : :
439 : : double P0[2], P1[2], AA[2], BB[2];
440 [ # # ][ # # ]: 0 : CubitVector absnorm(fabs(norm.x()), fabs(norm.y()), fabs(norm.z()));
[ # # ][ # # ]
441 [ # # ][ # # ]: 0 : if (absnorm.x() >= absnorm.y() && absnorm.x() >= absnorm.z())
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
442 : : {
443 [ # # ][ # # ]: 0 : P0[0] = pt0->coordinates().y(); P0[1] = pt0->coordinates().z();
[ # # ][ # # ]
444 [ # # ][ # # ]: 0 : P1[0] = pt1->coordinates().y(); P1[1] = pt1->coordinates().z();
[ # # ][ # # ]
445 [ # # ][ # # ]: 0 : AA[0] = aa.y(); AA[1] = aa.z();
446 [ # # ][ # # ]: 0 : BB[0] = bb.y(); BB[1] = bb.z();
447 : : }
448 [ # # ][ # # ]: 0 : else if (absnorm.y() >= absnorm.x() && absnorm.y() >= absnorm.z())
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
449 : : {
450 [ # # ][ # # ]: 0 : P0[0] = pt0->coordinates().z(); P0[1] = pt0->coordinates().x();
[ # # ][ # # ]
451 [ # # ][ # # ]: 0 : P1[0] = pt1->coordinates().z(); P1[1] = pt1->coordinates().x();
[ # # ][ # # ]
452 [ # # ][ # # ]: 0 : AA[0] = aa.z(); AA[1] = aa.x();
453 [ # # ][ # # ]: 0 : BB[0] = bb.z(); BB[1] = bb.x();
454 : : }
455 : : else
456 : : {
457 [ # # ][ # # ]: 0 : P0[0] = pt0->coordinates().x(); P0[1] = pt0->coordinates().y();
[ # # ][ # # ]
458 [ # # ][ # # ]: 0 : P1[0] = pt1->coordinates().x(); P1[1] = pt1->coordinates().y();
[ # # ][ # # ]
459 [ # # ][ # # ]: 0 : AA[0] = aa.x(); AA[1] = aa.y();
460 [ # # ][ # # ]: 0 : BB[0] = bb.x(); BB[1] = bb.y();
461 : : }
462 : :
463 : : double QQ[4], s;
464 [ # # ]: 0 : int ninter = intersect_2D_segments(P0, P1, AA, BB, QQ);
465 : :
466 [ # # ]: 0 : if (ninter != 1)
467 : : {
468 : 0 : does_intersect = CUBIT_FALSE;
469 : 0 : return CUBIT_SUCCESS;
470 : : }
471 : 0 : does_intersect = CUBIT_TRUE;
472 : :
473 : 0 : double dx = P1[0] - P0[0];
474 : 0 : double dy = P1[1] - P0[1];
475 [ # # ]: 0 : if (fabs(dx) > fabs(dy))
476 : 0 : s = (QQ[0] - P0[0]) / dx;
477 : : else
478 : 0 : s = (QQ[1] - P0[1]) / dy;
479 : :
480 [ # # ][ # # ]: 0 : qq = pt0->coordinates() + s * (pt1->coordinates() - pt0->coordinates());
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
481 : :
482 : 0 : return CUBIT_SUCCESS;
483 : : }
484 : :
485 : : //======================================================================
486 : : // Function: boundary_edge_points (PUBLIC)
487 : : // Description: return the oriented endpoints of a facet edge assuming
488 : : // the edge is at the frontier of the facets (only one adjacency)
489 : : // Author: sjowen
490 : : // Date: 8/00
491 : : //======================================================================
492 : 0 : void CubitFacetEdge::boundary_edge_points( CubitPoint * &pt0,
493 : : CubitPoint * &pt1,
494 : : int tool_id )
495 : : {
496 [ # # ]: 0 : if(num_adj_facets() == 0)
497 : : {
498 : 0 : pt0 = point(0);
499 : 0 : pt1 = point(1);
500 : 0 : return;
501 : : }
502 : 0 : CubitFacet *facet = NULL;
503 [ # # ]: 0 : if (tool_id != 0)
504 : : {
505 : : int ii;
506 : 0 : int found = 0;
507 [ # # ]: 0 : DLIList <CubitFacet *> adj_facets;
508 [ # # ]: 0 : facets(adj_facets);
509 [ # # ][ # # ]: 0 : for (ii=0; ii<adj_facets.size() && !found; ii++)
[ # # ][ # # ]
510 : : {
511 [ # # ]: 0 : facet = adj_facets.get_and_step();
512 [ # # ][ # # ]: 0 : if (facet->tool_id() == tool_id)
513 : 0 : found = 1;
514 : : }
515 [ # # ][ # # ]: 0 : assert(found);
516 : : }
517 : : else
518 : : {
519 : 0 : facet = adj_facet(0);
520 [ # # ]: 0 : if (!facet) facet = adj_facet(1);
521 [ # # ]: 0 : assert(facet != 0);
522 : : }
523 : :
524 : 0 : int index = facet->edge_index( this );
525 : 0 : int use = facet->edge_use( index );
526 [ # # ]: 0 : if (use == 1) {
527 : 0 : pt0 = point(0);
528 : 0 : pt1 = point(1);
529 : : }
530 : : else {
531 : 0 : pt0 = point(1);
532 : 0 : pt1 = point(0);
533 : : }
534 : : }
535 : :
536 : : //======================================================================
537 : : // Function: dist_to_edge
538 : : // Description: return the distance from the point to this edge
539 : : // Author: sjowen
540 : : // Date: 2/01
541 : : // Corrected by JFowler 5/03
542 : : //======================================================================
543 : 0 : double CubitFacetEdge::dist_to_edge(
544 : : const CubitVector &this_point,
545 : : CubitVector &close_point,
546 : : CubitBoolean &outside_edge )
547 : : {
548 : 0 : double dist = 0.0;
549 [ # # ][ # # ]: 0 : CubitVector p0 = point(0)->coordinates();
550 [ # # ][ # # ]: 0 : CubitVector p1 = point(1)->coordinates();
551 [ # # ]: 0 : CubitVector edge_vec( p1, p0 );
552 [ # # ]: 0 : CubitVector point_vec( this_point, p0 );
553 : : double edge_length;
554 [ # # ]: 0 : edge_length = edge_vec.normalize();
555 [ # # ]: 0 : double dist_on_edge = edge_vec % point_vec;
556 [ # # ]: 0 : if (dist_on_edge < 0.0e0)
557 : : {
558 [ # # ]: 0 : close_point = p0;
559 : 0 : outside_edge = CUBIT_TRUE;
560 : : }
561 [ # # ]: 0 : else if (dist_on_edge > edge_length)
562 : : {
563 [ # # ]: 0 : close_point = p1;
564 : 0 : outside_edge = CUBIT_TRUE;
565 : : }
566 : : else
567 : : {
568 [ # # ][ # # ]: 0 : close_point = p0 - edge_vec * dist_on_edge;
[ # # ]
569 : 0 : outside_edge = CUBIT_FALSE;
570 : : }
571 [ # # ]: 0 : dist = close_point.distance_between( this_point );
572 : 0 : return dist;
573 : : }
574 : :
575 : : /*//======================================================================
576 : : // Function: dist_to_edge
577 : : // Description: return the distance from the point to this edge
578 : : // Author: sjowen
579 : : // Date: 2/01
580 : : //======================================================================
581 : : double CubitFacetEdge::dist_to_edge(
582 : : CubitVector &this_point,
583 : : CubitVector &close_point,
584 : : CubitBoolean &outside_edge )
585 : : {
586 : : double dist = 0.0;
587 : : CubitVector p0 = point(0)->coordinates();
588 : : CubitVector p1 = point(1)->coordinates();
589 : : CubitVector edge_vec( p1, p0 );
590 : : CubitVector point_vec( this_point, p0 );
591 : : point_vec.normalize();
592 : : double dist_on_edge = edge_vec % point_vec;
593 : : double edge_length;
594 : : edge_length = edge_vec.normalize();
595 : : if (dist_on_edge < 0.0e0)
596 : : {
597 : : close_point = p0;
598 : : outside_edge = CUBIT_TRUE;
599 : : }
600 : : else if (dist_on_edge > edge_length)
601 : : {
602 : : close_point = p1;
603 : : outside_edge = CUBIT_TRUE;
604 : : }
605 : : else
606 : : {
607 : : close_point = p0 + edge_vec * dist_on_edge;
608 : : outside_edge = CUBIT_FALSE;
609 : : }
610 : : dist = close_point.distance_between( this_point );
611 : : return dist;
612 : : }
613 : : */
614 : : //======================================================================
615 : : // Function: proj_to_line
616 : : // Description: project the point to a line defined by the edge
617 : : // Author: sjowen
618 : : // Date: 2/01
619 : : //======================================================================
620 : 0 : CubitStatus CubitFacetEdge::proj_to_line(
621 : : const CubitVector &this_point,
622 : : CubitVector &proj_point )
623 : : {
624 : 0 : CubitStatus stat = CUBIT_SUCCESS;
625 [ # # ][ # # ]: 0 : CubitVector p0 = point(0)->coordinates();
626 [ # # ][ # # ]: 0 : CubitVector p1 = point(1)->coordinates();
627 [ # # ]: 0 : CubitVector edge_vec( p0,p1 );
628 [ # # ]: 0 : CubitVector point_vec( p0, this_point );
629 [ # # ]: 0 : edge_vec.normalize();
630 [ # # ]: 0 : double dist_on_edge = edge_vec % point_vec;
631 [ # # ][ # # ]: 0 : proj_point = p0 + (edge_vec * dist_on_edge);
[ # # ]
632 : :
633 : 0 : return stat;
634 : : }
635 : :
636 : : //======================================================================
637 : : // Function: edge_tangent
638 : : // Description: return tangent at a point on the edge
639 : : // Author: sjowen
640 : : // Date: 2/01
641 : : //======================================================================
642 : 330 : CubitStatus CubitFacetEdge::edge_tangent(
643 : : const CubitVector &/*point_on_edge*/,
644 : : CubitVector &tangent )
645 : : {
646 : 330 : CubitStatus stat = CUBIT_SUCCESS;
647 [ + - ][ + - ]: 330 : tangent = point(1)->coordinates() -
[ + - ]
648 [ + - ]: 660 : point(0)->coordinates();
649 : 330 : tangent.normalize();
650 : 330 : return stat;
651 : : }
652 : :
653 : : //======================================================================
654 : : // Function: edge_tangent
655 : : // Description: return curvature at a point on the edge
656 : : // Author: sjowen
657 : : // Date: 2/01
658 : : //======================================================================
659 : 0 : CubitStatus CubitFacetEdge::edge_curvature(
660 : : const CubitVector &/*point_on_edge*/,
661 : : CubitVector &curvature,
662 : : CubitFacetEdge *closest_edge )
663 : : {
664 [ # # ][ # # ]: 0 : CubitVector vec_ba, vec_ca, center_point;
[ # # ]
665 : :
666 : : //if point(0) is middle point
667 [ # # ][ # # ]: 0 : if( closest_edge->other_point( point(0) ) )
[ # # ]
668 : : {
669 [ # # ][ # # ]: 0 : center_point = point(0)->coordinates();
[ # # ]
670 [ # # ][ # # ]: 0 : vec_ba = closest_edge->point(0)->coordinates() - center_point;
[ # # ][ # # ]
671 [ # # ][ # # ]: 0 : vec_ca = point(1)->coordinates() - center_point;
[ # # ][ # # ]
672 : : }
673 : : //if point(1) is middle point
674 [ # # ][ # # ]: 0 : else if( closest_edge->other_point( point(1) ) )
[ # # ]
675 : : {
676 [ # # ][ # # ]: 0 : center_point = point(1)->coordinates();
[ # # ]
677 [ # # ][ # # ]: 0 : vec_ba = point(0)->coordinates() - center_point;
[ # # ][ # # ]
678 [ # # ][ # # ]: 0 : vec_ca = closest_edge->point(1)->coordinates() - center_point;
[ # # ][ # # ]
679 : : }
680 : : else
681 : 0 : assert(0);
682 : :
683 : : // Squares of lengths of the edges incident to `a'.
684 [ # # ]: 0 : double ba_length = vec_ba.length_squared();
685 [ # # ]: 0 : double ca_length = vec_ca.length_squared();
686 : :
687 : : // Cross product of these edges.
688 : : // (Take your chances with floating-point roundoff.)
689 [ # # ]: 0 : CubitVector cross_bc = vec_ba * vec_ca;
690 : :
691 : : // Calculate the denominator of the formulae.
692 [ # # ]: 0 : double temp_dbl = cross_bc % cross_bc;
693 [ # # ]: 0 : CubitVector circle_center(0.0,0.0,0.0);
694 [ # # ]: 0 : if(fabs(temp_dbl) > CUBIT_DBL_MIN){
695 : 0 : double denominator = 0.5 / (temp_dbl);
696 [ # # ]: 0 : assert(denominator != 0.0);
697 : :
698 : : // Calculate offset (from `a') of circumcenter.
699 [ # # ][ # # ]: 0 : circle_center = (ba_length * vec_ca - ca_length * vec_ba) * cross_bc;
[ # # ][ # # ]
[ # # ]
700 [ # # ]: 0 : circle_center *= denominator;
701 : :
702 : : //store radius
703 [ # # ]: 0 : double radius = circle_center.length();
704 [ # # ]: 0 : circle_center.normalize();
705 [ # # ]: 0 : circle_center /= radius;
706 : : }
707 [ # # ]: 0 : curvature = circle_center;
708 : :
709 : 0 : return CUBIT_SUCCESS;
710 : : }
711 : :
712 : : //======================================================================
713 : : // Function: length
714 : : // Description: return length of an edge
715 : : // Author: sjowen
716 : : // Date: 2/01
717 : : //======================================================================
718 : 946 : double CubitFacetEdge::length()
719 : : {
720 [ + - ][ + - ]: 946 : CubitVector start = point(0)->coordinates();
721 [ + - ][ + - ]: 946 : CubitVector end = point(1)->coordinates();
722 [ + - ]: 946 : return start.distance_between( end );
723 : : }
724 : :
725 : : //======================================================================
726 : : // Function: length
727 : : // Description: return length of an edge
728 : : // Author: jakraft
729 : : // Date: 3/04
730 : : //======================================================================
731 : 0 : CubitVector CubitFacetEdge::position_from_fraction( double f )
732 : : {
733 [ # # ][ # # ]: 0 : return (1.0 - f) * point(0)->coordinates() +
[ # # ]
734 [ # # ][ # # ]: 0 : f * point(1)->coordinates();
735 : : }
736 : :
737 : : //======================================================================
738 : : // Function: other_point
739 : : // Description: return the other point on the edge
740 : : // Author: sjowen
741 : : // Date: 4/01
742 : : //======================================================================
743 : 4532 : CubitPoint* CubitFacetEdge::other_point( CubitPoint *point_ptr )
744 : : {
745 [ + + ]: 4532 : if (point(0) == point_ptr)
746 : 3289 : return point(1);
747 [ + - ]: 1243 : if(point(1) == point_ptr)
748 : 1243 : return point(0);
749 : 0 : return NULL;
750 : : }
751 : :
752 : : //======================================================================
753 : : // Function: get_parents
754 : : // Description: return entities attached to this edge
755 : : // Author: sjowen
756 : : // Date: 4/01
757 : : //======================================================================
758 : 10318 : void CubitFacetEdge::get_parents(DLIList<FacetEntity *> &facet_list)
759 : : {
760 [ + - ]: 10318 : DLIList<CubitFacet *> cf_list;
761 [ + - ]: 10318 : facets( cf_list );
762 [ + - ][ + + ]: 30668 : for (int ii=0; ii<cf_list.size(); ii++)
763 [ + - ][ + - ]: 30668 : facet_list.append(cf_list.get_and_step());
[ + - ]
764 : 10318 : }
765 : :
766 : : //======================================================================
767 : : // Function: other_facet
768 : : // Description: return the other facet on the edge (assumes
769 : : // max two facets per edge)
770 : : // Author: sjowen
771 : : // Date: 5/01
772 : : //======================================================================
773 : 1936 : CubitFacet *CubitFacetEdge::other_facet( CubitFacet *facet_ptr )
774 : : {
775 [ + - ]: 1936 : DLIList<CubitFacet *> cf_list;
776 [ + - ]: 1936 : facets( cf_list );
777 [ + - ][ - + ]: 1936 : assert(cf_list.size() < 3);
778 : 1936 : CubitFacet *adj_facet = NULL;
779 [ + - ][ + - ]: 1936 : if (cf_list.size() > 0)
780 : : {
781 [ + - ]: 1936 : adj_facet = cf_list.get_and_step();
782 [ + + ]: 1936 : if (adj_facet == facet_ptr)
783 : : {
784 [ + - ][ + - ]: 968 : if (cf_list.size() == 2)
785 : : {
786 [ + - ]: 968 : adj_facet = cf_list.get();
787 : : }
788 : : }
789 : : }
790 [ + - ]: 1936 : return adj_facet;
791 : : }
792 : :
793 : : //======================================================================
794 : : // Function: other_facet_on_surf
795 : : // Description: return the other facet on the edge that has the
796 : : // same tool id. (finds the first occurrence
797 : : // of the same tool id of the adjacent facets
798 : : // to this edge that is not facet_ptr)
799 : : // Author: sjowen
800 : : // Date: 5/01
801 : : //======================================================================
802 : 2244 : CubitFacet *CubitFacetEdge::other_facet_on_surf( CubitFacet *facet_ptr )
803 : : {
804 [ - + ]: 2244 : assert(facet_ptr != 0);
805 [ + - ]: 2244 : int tool_id = facet_ptr->tool_id();
806 [ + - ]: 2244 : DLIList<CubitFacet *> cf_list;
807 [ + - ]: 2244 : facets( cf_list );
808 : 2244 : CubitFacet *adj_facet = NULL;
809 : 2244 : int found = 0;
810 [ + - ][ + + ]: 6160 : for (int ii=0; ii<cf_list.size() && !found; ii++)
[ + + ][ + + ]
811 : : {
812 [ + - ]: 3916 : adj_facet = cf_list.get_and_step();
813 [ + + ]: 3916 : if (adj_facet != facet_ptr)
814 : : {
815 [ + - ][ + + ]: 2024 : if (adj_facet->tool_id() == tool_id)
816 : 968 : found = 1;
817 : : }
818 : : }
819 [ + + ]: 2244 : if (!found)
820 : 1276 : adj_facet = NULL;
821 [ + - ]: 2244 : return adj_facet;
822 : : }
823 : :
824 : : //======================================================================
825 : : // Function: num_adj_facets_on_surf
826 : : // Description: count the number of adjacent facets to this edge that
827 : : // have the specified tool id
828 : : // Author: sjowen
829 : : // Date: 5/01
830 : : //======================================================================
831 : 858 : int CubitFacetEdge::num_adj_facets_on_surf( int tool_id )
832 : : {
833 : :
834 [ + - ]: 858 : DLIList<CubitFacet *> cf_list;
835 [ + - ]: 858 : facets( cf_list );
836 : 858 : CubitFacet *adj_facet = NULL;
837 : 858 : int nfacets = 0;
838 [ + - ][ + + ]: 2508 : for (int ii=0; ii<cf_list.size(); ii++)
839 : : {
840 [ + - ]: 1650 : adj_facet = cf_list.get_and_step();
841 [ + - ][ + + ]: 1650 : if (adj_facet->tool_id() == tool_id)
842 : 1386 : nfacets++;
843 : : }
844 [ + - ]: 858 : return nfacets;
845 : : }
846 : :
847 : : //======================================================================
848 : : // Function: adj_facet_on_surf
849 : : // Description: return the first facet on the adjacent facet list with
850 : : // the indicated tool id
851 : : // Author: sjowen
852 : : // Date: 5/01
853 : : //======================================================================
854 : 330 : CubitFacet *CubitFacetEdge::adj_facet_on_surf( int tool_id )
855 : : {
856 [ + - ]: 330 : DLIList<CubitFacet *> cf_list;
857 [ + - ]: 330 : facets( cf_list );
858 : 330 : CubitFacet *adj_facet = NULL;
859 : 330 : int found = 0;
860 [ + - ][ + + ]: 792 : for (int ii=0; ii<cf_list.size() && !found; ii++)
[ + + ][ + + ]
861 : : {
862 [ + - ]: 462 : adj_facet = cf_list.get_and_step();
863 [ + - ][ + + ]: 462 : if (adj_facet->tool_id() == tool_id)
864 : : {
865 : 330 : found = 1;
866 : : }
867 : : }
868 [ - + ]: 330 : if (!found)
869 : 0 : adj_facet = NULL;
870 [ + - ]: 330 : return adj_facet;
871 : : }
872 : :
873 : : //======================================================================
874 : : // Function: contains
875 : : // Description: determines if point is contained in edge
876 : : // Author: sjowen
877 : : // Date: 5/01
878 : : //======================================================================
879 : 6732 : CubitBoolean CubitFacetEdge::contains( CubitPoint *point_ptr )
880 : : {
881 [ + + ][ + + ]: 6732 : if (point(0) == point_ptr || point(1) == point_ptr)
[ + + ]
882 : 4488 : return CUBIT_TRUE;
883 : 2244 : return CUBIT_FALSE;
884 : : }
885 : :
886 : : //======================================================================
887 : : // Function: draw
888 : : // Description: draw the edge
889 : : // Author: sjowen
890 : : // Date: 5/01
891 : : //======================================================================
892 : 0 : void CubitFacetEdge::debug_draw(int color, int flush, int /*draw_uv*/)
893 : : {
894 [ # # ]: 0 : if ( color == -1 )
895 : 0 : color = CUBIT_RED_INDEX;
896 : 0 : GfxDebug::draw_facet_edge(this, color);
897 [ # # ]: 0 : GfxDebug::draw_point( point(0)->coordinates(), color );
898 [ # # ]: 0 : GfxDebug::draw_point( point(1)->coordinates(), color );
899 [ # # ]: 0 : if ( flush )
900 : 0 : GfxDebug::flush();
901 : 0 : }
902 : :
903 : : //======================================================================
904 : : // Function: shared_point (PUBLIC)
905 : : // Description: find the common point
906 : : // Author: sjowen
907 : : // Date: 10/02
908 : : //======================================================================
909 : 66 : CubitPoint *CubitFacetEdge::shared_point( CubitFacetEdge *edge_ptr )
910 : : {
911 : 66 : CubitPoint *pA = this->point(0);
912 : 66 : CubitPoint *pB = this->point(1);
913 : 66 : CubitPoint *pC = edge_ptr->point(0);
914 : 66 : CubitPoint *pD = edge_ptr->point(1);
915 : 66 : CubitPoint *pShared = NULL;
916 [ + - ][ + + ]: 66 : if (pA == pC || pA == pD)
917 : : {
918 : 55 : pShared = pA;
919 : : }
920 [ - + ][ # # ]: 11 : else if (pB == pC || pB == pD)
921 : : {
922 : 11 : pShared = pB;
923 : : }
924 : 66 : return pShared;
925 : : }
926 : :
927 : :
928 : : //======================================================================
929 : : // Function: bounding_box (PUBLIC)
930 : : // Description: return the bounding box of the edge
931 : : // Author: sjowen
932 : : // Date: 10/02
933 : : //======================================================================
934 : 0 : CubitBox CubitFacetEdge::bounding_box( )
935 : : {
936 [ # # ]: 0 : CubitPoint *p1 = point (0);
937 [ # # ]: 0 : CubitPoint *p2 = point (1);
938 : :
939 [ # # ][ # # ]: 0 : CubitVector bbox_min, bbox_max;
940 [ # # ][ # # ]: 0 : bbox_min.x(CUBIT_MIN(p1->x(),p2->x()));
[ # # ][ # # ]
[ # # ][ # # ]
941 [ # # ][ # # ]: 0 : bbox_min.y(CUBIT_MIN(p1->y(),p2->y()));
[ # # ][ # # ]
[ # # ][ # # ]
942 [ # # ][ # # ]: 0 : bbox_min.z(CUBIT_MIN(p1->z(),p2->z()));
[ # # ][ # # ]
[ # # ][ # # ]
943 [ # # ][ # # ]: 0 : bbox_max.x(CUBIT_MAX(p1->x(),p2->x()));
[ # # ][ # # ]
[ # # ][ # # ]
944 [ # # ][ # # ]: 0 : bbox_max.y(CUBIT_MAX(p1->y(),p2->y()));
[ # # ][ # # ]
[ # # ][ # # ]
945 [ # # ][ # # ]: 0 : bbox_max.z(CUBIT_MAX(p1->z(),p2->z()));
[ # # ][ # # ]
[ # # ][ # # ]
946 [ # # ]: 0 : CubitBox edge_box(bbox_min,bbox_max);
947 : 0 : return edge_box;
948 : : }
949 : :
950 : : //======================================================================
951 : : //
952 : : // The following 2 function really belong in IntersectionTool in the util
953 : : // directory. Because of the current distribution policy it was not easy
954 : : // to get it into there without breaking the build. These should be
955 : : // removed at a later time and integrated with the IntersectionTool
956 : : //
957 : : //=======================================================================
958 : : // Intersection of two line segments in 2D. Intersect [P0,P1] with [P2,P3]
959 : : // Return value is zero if there is no intersection, 1 if there is a unique
960 : : // intersection, and 2 if the 2 segments overlap and the intersection set is
961 : : // a segment itself. The return value is the number of valid entries*2 in
962 : : // the array qq.
963 : 0 : int CubitFacetEdge::intersect_2D_segments( double P0[2], double P1[2],
964 : : double P2[2], double P3[2],
965 : : double qq[4] )
966 : : {
967 : :
968 : : double D0[2], D1[2];
969 : 0 : D0[0] = P1[0] - P0[0]; D0[1] = P1[1] - P0[1];
970 : 0 : D1[0] = P3[0] - P2[0]; D1[1] = P3[1] - P2[1];
971 : :
972 : : // segments P0 + s * D0 for s in [0,1],
973 : : // P2 + t * D1 for t in [0,1]
974 : :
975 : 0 : double sqr_epsilon = DBL_EPSILON * DBL_EPSILON;
976 : : double E[2];
977 : 0 : E[0] = P2[0] - P0[0];
978 : 0 : E[1] = P2[1] - P0[1];
979 : :
980 : 0 : double kross = D0[0] * D1[1] - D0[1] * D1[0];
981 : 0 : double sqr_kross = kross * kross;
982 : 0 : double sqr_len0 = D0[0] * D0[0] + D0[1] * D0[1];
983 : 0 : double sqr_len1 = D1[0] * D1[0] + D1[1] * D1[1];
984 [ # # ]: 0 : if (sqr_kross > sqr_epsilon * sqr_len0 * sqr_len1)
985 : : {
986 : : // lines of the segment are not parallel
987 : :
988 : 0 : double s = (E[0] * D1[1] - E[1] * D1[0]) / kross;
989 [ # # ][ # # ]: 0 : if (s < 0.0 || s > 1.0)
990 : : {
991 : : // intersection of lines is not a point on segment P0 + s * D0
992 : 0 : return 0;
993 : : }
994 : :
995 : 0 : double t = (E[0] * D0[1] - E[1] * D0[0]) / kross;
996 [ # # ][ # # ]: 0 : if (t < 0.0 || t > 1.0)
997 : : {
998 : : // intersection of lines is not a point on segment P1 + t * D1
999 : 0 : return 0;
1000 : : }
1001 : :
1002 : : // intersection of lines is a point on each segment
1003 : :
1004 : 0 : qq[0] = P0[0] + s * D0[0];
1005 : 0 : qq[1] = P0[1] + s * D0[1];
1006 : 0 : return 1;
1007 : : }
1008 : :
1009 : : // lines of the segments are parallel
1010 : :
1011 : 0 : double sqr_lenE = E[0] * E[0] + E[1] * E[1];
1012 : 0 : kross = E[0] * D0[1] - E[1] * D0[1];
1013 : 0 : sqr_kross = kross * kross;
1014 [ # # ]: 0 : if (sqr_kross > sqr_epsilon * sqr_len0 * sqr_lenE)
1015 : : {
1016 : : // lines of the segments are different
1017 : 0 : return 0;
1018 : : }
1019 : :
1020 : : // lines of the segment are the same. Need to test for overlap of segments
1021 : :
1022 : 0 : double s0 = (D0[0] * E[0] + D0[1] * E[1]) / sqr_len0;
1023 : 0 : double s1 = s0 + (D0[0] * D1[0] + D0[1] * D1[1]) / sqr_len0;
1024 [ # # ]: 0 : double smin = CUBIT_MIN(s0, s1);
1025 [ # # ]: 0 : double smax = CUBIT_MAX(s0, s1);
1026 : :
1027 : : double w[2];
1028 [ # # ]: 0 : int imax = intersect_intervals(0.0, 1.0, smin, smax, w);
1029 [ # # ]: 0 : for (int i=0; i<imax; i++)
1030 : : {
1031 : 0 : qq[i*2] = P0[0] + w[i] * D0[0];
1032 : 0 : qq[i*2+1] = P0[1] + w[i] * D0[1];
1033 : : }
1034 : 0 : return imax;
1035 : : }
1036 : :
1037 : : // The intersection of two intervals [u0,u1] and [v0,v1], where u0<u1 and v0<v1.
1038 : : // Return value is zero if intervals do not intersect; 1 if they intersect at
1039 : : // a single point, in which case w[0] contains that point; or 2 if they intersect
1040 : : // in an interval whose endpoints are stored in w[0] and w[1]
1041 : 0 : int CubitFacetEdge::intersect_intervals( double u0, double u1,
1042 : : double v0, double v1,
1043 : : double w[2] )
1044 : : {
1045 [ # # ][ # # ]: 0 : if (u1 < v0 || u0 > v1)
1046 : 0 : return 0;
1047 : :
1048 [ # # ]: 0 : if (u1 > v0)
1049 : : {
1050 [ # # ]: 0 : if (u0 < v1)
1051 : : {
1052 [ # # ]: 0 : if (u0 < v0)
1053 : 0 : w[0] = v0;
1054 : : else
1055 : 0 : w[0] = u0;
1056 [ # # ]: 0 : if (u1 > v1)
1057 : 0 : w[1] = v1;
1058 : : else
1059 : 0 : w[1] = u1;
1060 : 0 : return 2;
1061 : : }
1062 : : else
1063 : : {
1064 : 0 : w[0] = u0;
1065 : 0 : return 1;
1066 : : }
1067 : : }
1068 : : else
1069 : : {
1070 : 0 : w[0] = u1;
1071 : 0 : return 1;
1072 : : }
1073 : : return 0;
1074 : : }
1075 : :
1076 : : //======================================================================
1077 : : // Function: add_facets (PUBLIC)
1078 : : // Description: add this edge to its adjacent facets
1079 : : // Author: sjowen
1080 : : // Date: 04/04
1081 : : //======================================================================
1082 : 0 : void CubitFacetEdge::add_facets( )
1083 : : {
1084 [ # # ]: 0 : CubitPoint *p0 = point (0);
1085 [ # # ]: 0 : CubitPoint *p1 = point (1);
1086 : :
1087 [ # # ]: 0 : DLIList<CubitFacet *> adj_facets;
1088 [ # # ]: 0 : p0->shared_facets(p1, adj_facets );
1089 : :
1090 : : CubitFacet *facet;
1091 : : int ii;
1092 [ # # ][ # # ]: 0 : for(ii=0; ii<adj_facets.size(); ii++)
1093 : : {
1094 [ # # ]: 0 : facet = adj_facets.get_and_step();
1095 [ # # ]: 0 : facet->add_edge( this );
1096 [ # # ]: 0 : }
1097 : 0 : }
1098 : :
1099 : 0 : double CubitFacetEdge::angle_between_facets()
1100 : : {
1101 [ # # ]: 0 : CubitFacet *facet0 = adj_facet(0);
1102 [ # # ]: 0 : CubitFacet *facet1 = adj_facet(1);
1103 : :
1104 : : // assumes facets are always oriented with outwards pointing normal
1105 : :
1106 [ # # ]: 0 : CubitVector n0 = facet0->normal();
1107 [ # # ]: 0 : CubitVector n1 = facet1->normal();
1108 : :
1109 : : // get orientation of edge with respect to facet0
1110 : :
1111 [ # # ]: 0 : CubitPoint *p0 = point(0);
1112 [ # # ]: 0 : CubitPoint *p1 = point(1);
1113 [ # # ]: 0 : CubitPoint *pnext = facet0->next_node(p0);
1114 [ # # ]: 0 : if (pnext != p1)
1115 : : {
1116 : 0 : pnext = p0;
1117 : 0 : p0 = p1;
1118 : 0 : p1 = pnext;
1119 : : }
1120 [ # # ][ # # ]: 0 : CubitVector evec = p1->coordinates() - p0->coordinates();
[ # # ]
1121 [ # # ]: 0 : evec.normalize();
1122 [ # # ][ # # ]: 0 : CubitVector cross = n0 * n1; cross.normalize();
1123 : : double angle;
1124 [ # # ]: 0 : double edot = evec % cross;
1125 [ # # ]: 0 : double ndot = n0 % n1;
1126 [ # # ]: 0 : if (ndot >= 1.0)
1127 : 0 : angle = 0.0;
1128 [ # # ]: 0 : else if (ndot <= -1.0)
1129 : 0 : angle = CUBIT_PI;
1130 : : else
1131 : 0 : angle = acos(ndot);
1132 [ # # ]: 0 : if (edot <= 0.0)
1133 : : {
1134 : 0 : angle = 2.0 * CUBIT_PI - angle;
1135 : : }
1136 : 0 : return angle;
1137 : : }
1138 : :
1139 : : // order edges in list beginning at start_point
1140 : : // report the endpoint
1141 : : // return CUBIT_SUCCESS if all edges are connected and ordered successfully
1142 : : // otherwise return CUBIT_FAILURE, in which case no changes are made
1143 : 66 : CubitStatus CubitFacetEdge::order_edge_list(DLIList<CubitFacetEdge*> &edge_list,
1144 : : CubitPoint *start_point,
1145 : : CubitPoint *&end_point)
1146 : : {
1147 : : int i;
1148 [ - + ]: 66 : assert(start_point);
1149 : :
1150 : 66 : end_point = NULL;
1151 : :
1152 : : // invalid input
1153 [ + - ][ - + ]: 66 : if (0 == edge_list.size())
1154 : 0 : return CUBIT_FAILURE;
1155 : :
1156 : : // simple case of a single edge - endpoitn
1157 [ + - ][ - + ]: 66 : if (1 == edge_list.size())
1158 : : {
1159 [ # # ][ # # ]: 0 : end_point = edge_list.get()->other_point(start_point);
1160 [ # # ]: 0 : return end_point ? CUBIT_SUCCESS : CUBIT_FAILURE;
1161 : : }
1162 : :
1163 [ + - ]: 66 : edge_list.reset();
1164 : :
1165 : : // note that a periodic/closed curve will fail
1166 : : // we could handle that case here if needed, but we may need more information
1167 : : // to know where to start and end the curve
1168 [ - + ]: 66 : if (NULL == start_point)
1169 : 0 : return CUBIT_FAILURE;
1170 : :
1171 : : // put edges in a set for faster searching
1172 [ + - ]: 66 : std::set<CubitFacetEdge *> edge_set;
1173 [ + - ][ + + ]: 286 : for (i=0; i<edge_list.size(); i++)
1174 [ + - ][ + - ]: 220 : edge_set.insert(dynamic_cast<CubitFacetEdge*> (edge_list.step_and_get()));
1175 : :
1176 : : // a vector for the ordered list
1177 [ + - ][ + - ]: 132 : std::vector<CubitFacetEdge*> ordered_edges;
1178 : :
1179 : : // find connected edges from the start point
1180 : 66 : CubitPoint *cur_pt = start_point;
1181 [ + + ]: 220 : do
1182 : : {
1183 : : // get edges connected to the current point and find the next edge
1184 [ + - ]: 220 : DLIList<CubitFacetEdge *> pt_edges;
1185 [ + - ]: 220 : cur_pt->edges(pt_edges);
1186 : :
1187 [ + - ]: 220 : std::set<CubitFacetEdge *>::iterator iter_found;
1188 : 220 : CubitFacetEdge *cur_edge = NULL;
1189 [ + - ][ + + ]: 561 : for (i=0; i<pt_edges.size() && !cur_edge; i++)
[ + + ][ + + ]
1190 : : {
1191 [ + - ]: 341 : CubitFacetEdge *tmp_edge = pt_edges.get_and_step();
1192 [ + - ]: 341 : iter_found = edge_set.find(tmp_edge);
1193 [ + - ][ + - ]: 341 : if ( iter_found != edge_set.end() )
[ + + ]
1194 : 220 : cur_edge = tmp_edge;
1195 : : }
1196 : :
1197 : : // if we don't find a connection before we empty the set
1198 : : // then not all the edges are connected -- return failure
1199 [ - + ]: 220 : if (NULL == cur_edge)
1200 : 0 : return CUBIT_FAILURE;
1201 : :
1202 : : // add the edge to the ordered list
1203 [ + - ]: 220 : ordered_edges.push_back( cur_edge );
1204 [ + - ]: 220 : edge_set.erase(iter_found);
1205 : :
1206 [ + - ][ + - ]: 220 : cur_pt = cur_edge->other_point(cur_pt);
[ + - ]
1207 : : }
1208 [ + - ]: 220 : while ( edge_set.size());
1209 : :
1210 [ + - ][ + - ]: 66 : if (ordered_edges.size() != (size_t)edge_list.size())
[ - + ]
1211 : 0 : return CUBIT_FAILURE;
1212 : :
1213 : : // store the edges in the correct order
1214 [ + - ]: 66 : edge_list.clean_out();
1215 : :
1216 [ + - ]: 66 : std::vector<CubitFacetEdge*>::iterator iter;
1217 [ + - ][ + - ]: 286 : for (iter=ordered_edges.begin(); iter!=ordered_edges.end(); iter++)
[ + - ][ + - ]
[ + + ]
1218 [ + - ][ + - ]: 220 : edge_list.append(*iter);
1219 : :
1220 : : // get the end point
1221 [ + - ][ + - ]: 66 : CubitFacetEdge *edge1 = edge_list[edge_list.size() - 1];
1222 [ + - ][ + - ]: 66 : CubitFacetEdge *edge2 = edge_list[edge_list.size() - 2];
1223 : :
1224 [ + - ][ + - ]: 66 : end_point = edge1->other_point( edge1->shared_point(edge2) );
1225 : :
1226 [ + - ]: 132 : return CUBIT_SUCCESS;
1227 : : }
1228 : :
1229 : 0 : CubitPoint *CubitFacetEdge::find_start_point_for_edge_list(DLIList<CubitFacetEdge*> edge_list)
1230 : : {
1231 : : // look for an edge with a point only connected to the one edge in the set
1232 : : //
1233 : : // TODO - this algorithm could be made more efficient by only checking one endpoint
1234 : : // per edge. The current implementation should match the order points were
1235 : : // checked in previous functionality.
1236 : : // If speed becomes an issue it could be reimplemented
1237 : : //
1238 : :
1239 [ # # ]: 0 : if (1 == edge_list.size())
1240 : : {
1241 : 0 : return edge_list[0]->point(0);
1242 : : }
1243 : :
1244 : : int i;
1245 : 0 : CubitPoint *start_point = NULL;
1246 [ # # ][ # # ]: 0 : for (i=0; i<edge_list.size() && (start_point == NULL); i++)
[ # # ]
1247 : : {
1248 [ # # ]: 0 : CubitFacetEdge *tmp_edge = edge_list.get_and_step();
1249 [ # # ]: 0 : DLIList<CubitFacetEdge*> pt_edges;
1250 [ # # ][ # # ]: 0 : tmp_edge->point(0)->edges(pt_edges);
1251 : :
1252 [ # # ]: 0 : pt_edges.intersect_unordered(edge_list);
1253 [ # # ][ # # ]: 0 : if (pt_edges.size() == 1)
1254 : : {
1255 [ # # ]: 0 : start_point = tmp_edge->point(0);
1256 : : }
1257 : : else
1258 : : {
1259 [ # # ]: 0 : pt_edges.clean_out();
1260 [ # # ][ # # ]: 0 : tmp_edge->point(1)->edges(pt_edges);
1261 [ # # ]: 0 : pt_edges.intersect_unordered(edge_list);
1262 [ # # ][ # # ]: 0 : if (pt_edges.size() == 1)
1263 : : {
1264 [ # # ]: 0 : start_point = tmp_edge->point(1);
1265 : : }
1266 : : }
1267 [ # # ]: 0 : }
1268 : 0 : return start_point;
1269 [ + - ][ + - ]: 6540 : }
1270 : :
1271 : : //EOF
|
