Branch data Line data Source code
1 : : //- Class: FacetDataUtil
2 : : //- Description: static library of general functions for querying and/or
3 : : //- modifying facet entities
4 : : //- Owner: Steve Owen
5 : : //- Checked by:
6 : : //- Version:
7 : :
8 : : #include "FacetDataUtil.hpp"
9 : : #include "FacetEvalTool.hpp"
10 : : #include "CubitFacet.hpp"
11 : : #include "CubitFacetEdge.hpp"
12 : : #include "CubitPoint.hpp"
13 : : #include "CubitPointData.hpp"
14 : : #include "CubitFacetData.hpp"
15 : : #include "CubitFacetEdgeData.hpp"
16 : : #include "CubitQuadFacet.hpp"
17 : : #include "GeometryDefines.h"
18 : : #include "ChollaDebug.hpp"
19 : : #include "GfxDebug.hpp"
20 : : #include "KDDTree.hpp"
21 : :
22 : : //===========================================================================
23 : : //Function Name: edges_by_count
24 : : //Description: find edges that are shared by "count" number of facets
25 : : //Author: bwhanks
26 : : //Date: 2003
27 : : //===========================================================================
28 : 0 : void FacetDataUtil::edges_by_count(DLIList<CubitFacet*> &facets,
29 : : unsigned int count,
30 : : DLIList<CubitFacetEdge*> &edges)
31 : : {
32 : : int i;
33 : : int j;
34 : :
35 : : // get a list of all edges from the facets
36 [ # # ]: 0 : DLIList<CubitFacetEdge*> edge_list;
37 [ # # ]: 0 : facets.reset();
38 [ # # ][ # # ]: 0 : for (i=facets.size(); i>0; i--)
39 : : {
40 [ # # ]: 0 : CubitFacet* p_facet = facets.get_and_step();
41 [ # # ]: 0 : for (j=0; j<3; j++)
42 : : {
43 [ # # ][ # # ]: 0 : edge_list.append(p_facet->edge(j));
44 : : }
45 : : }
46 : :
47 : : // reset edge marks
48 [ # # ]: 0 : edge_list.reset();
49 [ # # ][ # # ]: 0 : for (i=edge_list.size(); i>0; i--)
50 [ # # ][ # # ]: 0 : edge_list.get_and_step()->marked(0);
51 : :
52 : : // mark with the hit count
53 [ # # ]: 0 : edge_list.reset();
54 [ # # ][ # # ]: 0 : for (i=edge_list.size(); i>0; i--)
55 : : {
56 [ # # ]: 0 : CubitFacetEdge* p_edge = edge_list.get_and_step();
57 [ # # ][ # # ]: 0 : p_edge->marked(p_edge->marked() + 1);
58 : : }
59 : :
60 : : // create the output list of edges hit the number of times passed in
61 [ # # ]: 0 : edge_list.reset();
62 [ # # ][ # # ]: 0 : for (i=edge_list.size(); i>0; i--)
63 : : {
64 [ # # ]: 0 : CubitFacetEdge* p_edge = edge_list.get_and_step();
65 [ # # ][ # # ]: 0 : if (static_cast<unsigned>(p_edge->marked()) == count)
66 [ # # ]: 0 : edges.append(p_edge);
67 : : }
68 : :
69 : : // reset edge marks
70 [ # # ]: 0 : edge_list.reset();
71 [ # # ][ # # ]: 0 : for (i=edge_list.size(); i>0; i--)
72 [ # # ][ # # ]: 0 : edge_list.get_and_step()->marked(0);
[ # # ]
73 : 0 : }
74 : :
75 : : //===========================================================================
76 : : //Function Name: ordered_point_edge_bdry
77 : : //Description: get an ordered list of points and edges around the boundary
78 : : // of a set of facets. The output "chain" consists of point -
79 : : // edge - point - edge - ... around the boundary of the input
80 : : // facets. Having the ordered points as well as the edges
81 : : // helps provides sense information for the bounding edges.
82 : : //Author: bwhanks
83 : : //Date: 2003
84 : : //===========================================================================
85 : 0 : void FacetDataUtil::ordered_point_edge_bdry(DLIList<CubitFacet*> &facets,
86 : : DLIList<FacetEntity*> &point_edge_chain)
87 : : {
88 [ # # ][ # # ]: 0 : assert(point_edge_chain.size() == 0);
89 [ # # ]: 0 : DLIList<CubitFacetEdge*> unordered_edges;
90 : :
91 : 0 : int use_count = 1; // get boundary edges (used only once by facets in the facet list)
92 [ # # ]: 0 : FacetDataUtil::edges_by_count(facets, use_count, unordered_edges);
93 : :
94 : : // Adds the start_edge to the list of boundary_edges and then attempts to find the
95 : : // next edge by getting all of the edges belonging to the second point on this
96 : : // edge and searching for an unmarked (-1) boundary edge. Does this until it can't
97 : : // find another edge.
98 : :
99 : :
100 : : // mark all edges connected to points of boundary edges with -2
101 : : // NOTE: using negative marks, since as boundary edges are added to the ordered list they
102 : : // get marked with the order in which they are found
103 : :
104 : : // TODO - this seems very inefficient
105 : : int i;
106 : : int j;
107 : : int k;
108 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge*> pt_edge_list;
109 : : CubitFacetEdge* cur_edge;
110 [ # # ]: 0 : unordered_edges.reset();
111 [ # # ][ # # ]: 0 : for (i=unordered_edges.size(); i>0; i--)
112 : : {
113 [ # # ]: 0 : cur_edge = unordered_edges.get_and_step();
114 [ # # ]: 0 : for (j=0; j<2; j++)
115 : : {
116 [ # # ]: 0 : pt_edge_list.clean_out();
117 [ # # ][ # # ]: 0 : cur_edge->point(j)->edges(pt_edge_list);
118 [ # # ]: 0 : pt_edge_list.reset();
119 [ # # ][ # # ]: 0 : for (k=pt_edge_list.size(); k>0; k--)
120 [ # # ][ # # ]: 0 : pt_edge_list.get_and_step()->marked(-2);
121 : : }
122 : : }
123 : :
124 : : // mark all boundary edges with -1
125 [ # # ]: 0 : unordered_edges.reset();
126 [ # # ][ # # ]: 0 : for (i=unordered_edges.size(); i>0; i--)
127 [ # # ][ # # ]: 0 : unordered_edges.get_and_step()->marked(-1);
128 : :
129 : :
130 : :
131 : :
132 : : CubitPoint *edge_pt1, *edge_pt2, *pt2, *originalpt1;
133 : : CubitFacetEdge *start_edge;
134 : : CubitFacetEdge *this_edge;
135 : : CubitBoolean keepgoing;
136 : :
137 : 0 : int i_found_some = 0;
138 [ # # ]: 0 : unordered_edges.reset();
139 [ # # ]: 0 : start_edge = unordered_edges.get();
140 : :
141 : : // find the orientation of the first edge with respect to one of the facets
142 : : // then order the edges accordingly
143 [ # # ]: 0 : facets.reset();
144 : : int e_index;
145 : 0 : int start_edge_sense = 0;
146 [ # # ][ # # ]: 0 : for (i=facets.size(); i>0; i--)
147 : : {
148 [ # # ]: 0 : CubitFacet* p_facet = facets.get_and_step();
149 [ # # ][ # # ]: 0 : if ( (e_index = p_facet->edge_index(start_edge)) >= 0 )
150 : : {
151 [ # # ]: 0 : start_edge_sense = p_facet->edge_use(e_index);
152 : 0 : break;
153 : : }
154 : : }
155 : :
156 : :
157 [ # # ]: 0 : start_edge->marked(i_found_some);
158 : 0 : i_found_some += 1;
159 [ # # ]: 0 : if (1 == start_edge_sense)
160 : : {
161 [ # # ]: 0 : originalpt1 = start_edge->point(0);
162 [ # # ]: 0 : point_edge_chain.append(originalpt1);
163 [ # # ]: 0 : point_edge_chain.append(start_edge);
164 [ # # ]: 0 : pt2 = start_edge->point(1);
165 : : }
166 : : else
167 : : {
168 [ # # ]: 0 : assert(-1 == start_edge_sense);
169 [ # # ]: 0 : originalpt1 = start_edge->point(1);
170 [ # # ]: 0 : point_edge_chain.append(originalpt1);
171 [ # # ]: 0 : point_edge_chain.append(start_edge);
172 [ # # ]: 0 : pt2 = start_edge->point(0);
173 : : }
174 : :
175 : :
176 : : // Look for an edge having pt2 as a point and also being on the boundary.
177 [ # # ]: 0 : pt_edge_list.clean_out();
178 [ # # ]: 0 : pt2->edges(pt_edge_list);
179 : :
180 : 0 : keepgoing = CUBIT_TRUE;
181 [ # # ]: 0 : pt_edge_list.reset();
182 [ # # ]: 0 : while ( keepgoing == CUBIT_TRUE ) {
183 : 0 : keepgoing = CUBIT_FALSE;
184 [ # # ][ # # ]: 0 : for ( i = pt_edge_list.size(); i > 0; i-- ) {
185 [ # # ]: 0 : this_edge = pt_edge_list.get_and_step();
186 [ # # ][ # # ]: 0 : if ( this_edge->marked() == -1 ) {
187 : 0 : i_found_some += 1;
188 [ # # ]: 0 : this_edge->marked(i_found_some);
189 [ # # ]: 0 : point_edge_chain.append(pt2);
190 [ # # ]: 0 : point_edge_chain.append(this_edge);
191 : :
192 [ # # ]: 0 : edge_pt1 = this_edge->point(0);
193 [ # # ]: 0 : edge_pt2 = this_edge->point(1);
194 [ # # ]: 0 : if (pt2 == edge_pt1)
195 : : {
196 : 0 : pt2 = edge_pt2;
197 : : }
198 : : else
199 : : {
200 [ # # ]: 0 : assert(pt2 == edge_pt2);
201 : 0 : pt2 = edge_pt1;
202 : : }
203 : :
204 : :
205 : 0 : keepgoing = CUBIT_TRUE;
206 [ # # ]: 0 : pt_edge_list.clean_out();
207 [ # # ]: 0 : pt2->edges(pt_edge_list);
208 : 0 : break;
209 : : }
210 : : }
211 : : }
212 : :
213 : : // clear marks
214 [ # # ][ # # ]: 0 : for (i=unordered_edges.size(); i>0; i--)
215 : : {
216 [ # # ]: 0 : cur_edge = unordered_edges.get_and_step();
217 [ # # ]: 0 : for (j=0; j<2; j++)
218 : : {
219 [ # # ]: 0 : pt_edge_list.clean_out();
220 [ # # ][ # # ]: 0 : cur_edge->point(j)->edges(pt_edge_list);
221 [ # # ]: 0 : pt_edge_list.reset();
222 [ # # ][ # # ]: 0 : for (k=pt_edge_list.size(); k>0; k--)
223 [ # # ][ # # ]: 0 : pt_edge_list.get_and_step()->marked(0);
224 : : }
225 : : }
226 : :
227 [ # # ][ # # ]: 0 : assert(pt2 == originalpt1);
228 : 0 : }
229 : :
230 : : //===========================================================================
231 : : //Function Name: partial_chain
232 : : //Description: given an ordered point - edge boundary and a start and end
233 : : // point on the boundary, return the point - edge chain between
234 : : // the two points, inclusive.
235 : : //Author: bwhanks
236 : : //Date: 2003
237 : : //===========================================================================
238 : 0 : CubitStatus FacetDataUtil::partial_chain(DLIList<FacetEntity*> &point_edge_chain,
239 : : FacetEntity* point1,
240 : : FacetEntity* point2,
241 : : DLIList<FacetEntity*> &chain_between)
242 : : {
243 : : // make sure the edges_between list is empty
244 [ # # ]: 0 : assert(chain_between.size() == 0);
245 [ # # ]: 0 : assert(point1 != point2);
246 : :
247 : : // find the start point in the list
248 : 0 : int pt1_index = point_edge_chain.where_is_item(point1);
249 : :
250 [ # # ]: 0 : if (-1 == pt1_index)
251 : 0 : return CUBIT_FAILURE;
252 : :
253 : 0 : point_edge_chain.reset();
254 : 0 : point_edge_chain.step(pt1_index);
255 : :
256 : 0 : CubitBoolean b_done = CUBIT_FALSE;
257 [ # # ]: 0 : while (!b_done)
258 : : {
259 [ # # ]: 0 : if (point_edge_chain.get() == point2)
260 : 0 : b_done = CUBIT_TRUE;
261 : :
262 : 0 : chain_between.append(point_edge_chain.get_and_step());
263 : : }
264 : :
265 : 0 : return CUBIT_SUCCESS;
266 : : }
267 : :
268 : : //===========================================================================
269 : : //Function Name: get_facet_points
270 : : //Description: return a unique list of the points from a given set of facets
271 : : //Author: bwhanks
272 : : //Date: 2003
273 : : //===========================================================================
274 : 0 : void FacetDataUtil::get_facet_points(DLIList<CubitFacet*> &cubit_facets,
275 : : DLIList<CubitPoint*> &facet_points)
276 : : {
277 : : int i;
278 [ # # ][ # # ]: 0 : DLIList<CubitPoint*> cubit_points(cubit_facets.size() * 3);
279 [ # # ][ # # ]: 0 : for( i = cubit_facets.size(); i--; )
280 : : {
281 [ # # ]: 0 : CubitFacet* facet = cubit_facets.step_and_get();
282 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
283 : : {
284 [ # # ]: 0 : CubitPoint* pt = dynamic_cast<CubitPoint*>(facet->point(j));
285 [ # # ]: 0 : assert(!!pt);
286 [ # # ]: 0 : pt->marked(0);
287 [ # # ]: 0 : cubit_points.append(pt);
288 : : }
289 : : }
290 [ # # ][ # # ]: 0 : for( i = cubit_points.size(); i--; )
291 : : {
292 [ # # ]: 0 : CubitPoint* pt = cubit_points.step_and_get();
293 [ # # ][ # # ]: 0 : pt->marked( pt->marked() + 1);
294 : : }
295 [ # # ][ # # ]: 0 : for( i = cubit_points.size(); i--; )
296 : : {
297 [ # # ]: 0 : CubitPoint* pt = cubit_points.step_and_get();
298 [ # # ][ # # ]: 0 : pt->marked( pt->marked() - 1 );
299 [ # # ][ # # ]: 0 : if( pt->marked() > 0 )
300 [ # # ]: 0 : cubit_points.change_to(0);
301 : : }
302 [ # # ]: 0 : cubit_points.remove_all_with_value(0);
303 : :
304 [ # # ][ # # ]: 0 : facet_points = cubit_points;
305 : 0 : }
306 : :
307 : : //===========================================================================
308 : : //Function Name: get_boundary_points
309 : : //Description: return the boundary points from a list of facets (unordered)
310 : : // assumes edges exist on the facets
311 : : //Author: bwhanks
312 : : //Date: 2003
313 : : //===========================================================================
314 : 0 : void FacetDataUtil::get_boundary_points(DLIList<CubitFacet*> &facet_list,
315 : : DLIList<CubitPoint*> &point_list)
316 : : {
317 : : int ii, jj;
318 : : CubitPoint *pt;
319 : : CubitFacetEdge *edge;
320 : : CubitFacet *facet;
321 [ # # ]: 0 : DLIList<CubitFacetEdge *>edge_list;
322 [ # # ][ # # ]: 0 : for(ii=0; ii<facet_list.size(); ii++)
323 : : {
324 [ # # ]: 0 : facet = facet_list.get_and_step();
325 [ # # ]: 0 : for(jj=0; jj<3; jj++)
326 : : {
327 [ # # ]: 0 : edge = facet->edge( jj );
328 [ # # ][ # # ]: 0 : if (1 == edge->num_adj_facets())
329 : : {
330 [ # # ]: 0 : edge_list.append(edge);
331 [ # # ][ # # ]: 0 : edge->point( 0 )->marked( 0 );
332 [ # # ][ # # ]: 0 : edge->point( 1 )->marked( 0 );
333 : : }
334 : : }
335 : : }
336 [ # # ][ # # ]: 0 : for(ii=0; ii<edge_list.size(); ii++)
337 : : {
338 [ # # ]: 0 : edge = edge_list.get_and_step();
339 [ # # ]: 0 : pt = edge->point( 0 );
340 [ # # ][ # # ]: 0 : if (pt->marked() == 0)
341 : : {
342 [ # # ]: 0 : pt->marked( 1 );
343 [ # # ]: 0 : point_list.append( pt );
344 : : }
345 [ # # ]: 0 : pt = edge->point( 1 );
346 [ # # ][ # # ]: 0 : if (pt->marked() == 0)
347 : : {
348 [ # # ]: 0 : pt->marked( 1 );
349 [ # # ]: 0 : point_list.append( pt );
350 : : }
351 : : }
352 [ # # ][ # # ]: 0 : for(ii=0; ii<point_list.size(); ii++)
353 [ # # ][ # # ]: 0 : point_list.get_and_step()->marked(0);
[ # # ]
354 : :
355 : 0 : }
356 : :
357 : : //===========================================================================
358 : : //Function Name: get_boundary_edges
359 : : //Description: return an unordered list of edges at the boundary of a
360 : : // set of facets
361 : : //Author: sjowen
362 : : //Date: 1/19/2004
363 : : //===========================================================================
364 : 0 : void FacetDataUtil::get_boundary_edges(DLIList<CubitFacet*> &facet_list,
365 : : DLIList<CubitFacetEdge*> &edge_list)
366 : : {
367 : : int ii, jj;
368 : : CubitFacetEdge *edge;
369 : : CubitFacet *facet;
370 [ # # ][ # # ]: 0 : for(ii=0; ii<facet_list.size(); ii++)
371 : : {
372 [ # # ]: 0 : facet = facet_list.get_and_step();
373 [ # # ]: 0 : for(jj=0; jj<3; jj++)
374 : : {
375 [ # # ]: 0 : edge = facet->edge( jj );
376 [ # # ][ # # ]: 0 : if (1 == edge->num_adj_facets())
377 : : {
378 [ # # ]: 0 : edge_list.append(edge);
379 : : }
380 : : }
381 : : }
382 : 0 : }
383 : :
384 : : //===========================================================================
385 : : //Function Name: get_points
386 : : //Description: get a unique set of points from facets
387 : : //Author: sjowen
388 : : //Date: 9/11/03
389 : : //===========================================================================
390 : 44 : void FacetDataUtil::get_points(DLIList<CubitFacet*> &facet_list,
391 : : DLIList<CubitPoint*> &point_list)
392 : : {
393 : : CubitFacet *facet;
394 : : CubitPoint *pt;
395 : : int ii, jj;
396 [ + - ][ + + ]: 396 : for(ii=0; ii<facet_list.size(); ii++)
397 : : {
398 [ + - ]: 352 : facet = facet_list.get_and_step();
399 [ + - ][ + - ]: 352 : facet->point( 0 )->marked( 0 );
400 [ + - ][ + - ]: 352 : facet->point( 1 )->marked( 0 );
401 [ + - ][ + - ]: 352 : facet->point( 2 )->marked( 0 );
402 : : }
403 : :
404 [ + - ][ + + ]: 396 : for (ii=0; ii<facet_list.size(); ii++)
405 : : {
406 [ + - ]: 352 : facet = facet_list.get_and_step();
407 [ + + ]: 1408 : for(jj=0; jj<3; jj++)
408 : : {
409 [ + - ]: 1056 : pt = facet->point(jj);
410 [ + - ][ + + ]: 1056 : if (pt->marked() == 0)
411 : : {
412 [ + - ]: 352 : pt->marked(1);
413 [ + - ]: 352 : point_list.append(pt);
414 : : }
415 : : }
416 : : }
417 : :
418 [ + - ][ + + ]: 396 : for(ii=0; ii<facet_list.size(); ii++)
419 : : {
420 [ + - ]: 352 : facet = facet_list.get_and_step();
421 [ + - ][ + - ]: 352 : facet->point( 0 )->marked( 0 );
422 [ + - ][ + - ]: 352 : facet->point( 1 )->marked( 0 );
423 [ + - ][ + - ]: 352 : facet->point( 2 )->marked( 0 );
424 : : }
425 : 44 : }
426 : :
427 : : //===========================================================================
428 : : //Function Name: copy_facets
429 : : //Description: make a complete copy of facets, points and edges
430 : : //Note: copies edges and points with their "feature" flag
431 : : //Author: sjowen
432 : : //Date: 9/11/03
433 : : //===========================================================================
434 : 0 : void FacetDataUtil::copy_facets(DLIList<CubitFacet*> &old_facet_list,
435 : : DLIList<CubitFacet*> &new_facet_list,
436 : : DLIList<CubitPoint*> &new_point_list,
437 : : DLIList<CubitFacetEdge*> &new_edge_list)
438 : : {
439 : : // get a unique set of points from the facets
440 : :
441 [ # # ]: 0 : DLIList<CubitPoint *> old_point_list;
442 [ # # ]: 0 : get_points(old_facet_list, old_point_list);
443 [ # # ][ # # ]: 0 : CubitPoint **point_array = new CubitPoint* [old_point_list.size()];
[ # # ]
444 : :
445 : : //- copy the points
446 : :
447 : : int ii;
448 [ # # ]: 0 : old_point_list.reset();
449 : : CubitPoint *new_point, *the_point;
450 [ # # ][ # # ]: 0 : for(ii=0; ii<old_point_list.size(); ii++)
451 : : {
452 [ # # ]: 0 : the_point = old_point_list.get_and_step();
453 [ # # ][ # # ]: 0 : new_point = new CubitPointData( the_point->coordinates() );
[ # # ]
454 [ # # ]: 0 : the_point->marked( ii );
455 [ # # ]: 0 : new_point_list.append( new_point );
456 : 0 : point_array[ii] = new_point;
457 [ # # ][ # # ]: 0 : if (the_point->is_feature())
458 [ # # ]: 0 : new_point->set_as_feature();
459 : : }
460 : :
461 : : //- copy the facets
462 : :
463 : : int jj, idx;
464 : : CubitFacet *new_facet, *the_facet;
465 : : CubitPoint *points[3];
466 : :
467 [ # # ]: 0 : old_facet_list.reset();
468 [ # # ][ # # ]: 0 : for (ii=0; ii<old_facet_list.size(); ii++)
469 : : {
470 [ # # ]: 0 : the_facet = old_facet_list.get_and_step();
471 [ # # ]: 0 : for (jj=0; jj<3; jj++)
472 : : {
473 [ # # ][ # # ]: 0 : idx = the_facet->point(jj)->marked();
474 : 0 : points[jj] = point_array[idx];
475 : : }
476 [ # # ][ # # ]: 0 : new_facet = new CubitFacetData( points[0], points[1], points[2] );
477 [ # # ]: 0 : new_facet_list.append( new_facet );
478 : : }
479 : :
480 : : //- copy the edges
481 : :
482 : : int idx0, idx1;
483 : : CubitFacetEdge *new_edge;
484 : : CubitFacetEdge *old_edge;
485 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *>old_edge_list;
486 [ # # ]: 0 : get_edges(old_facet_list, old_edge_list);
487 [ # # ][ # # ]: 0 : for(ii=0; ii<old_edge_list.size(); ii++)
488 : : {
489 [ # # ]: 0 : old_edge = old_edge_list.get_and_step();
490 [ # # ][ # # ]: 0 : idx0 = old_edge->point(0)->marked();
491 [ # # ][ # # ]: 0 : idx1 = old_edge->point(1)->marked();
492 [ # # ][ # # ]: 0 : new_edge = new CubitFacetEdgeData( point_array[idx0], point_array[idx1] );
493 [ # # ][ # # ]: 0 : if (old_edge->is_feature())
494 [ # # ]: 0 : new_edge->set_as_feature();
495 [ # # ]: 0 : new_edge_list.append( new_edge );
496 : : }
497 : :
498 [ # # ][ # # ]: 0 : delete [] point_array;
499 : :
500 : 0 : }
501 : :
502 : : //===========================================================================
503 : : //Function Name: get_edges
504 : : //Description: Populates the edge list from the list of facets - creates the
505 : : // edges if necessary
506 : : //Author: sjowen
507 : : //Date: 9/11/03
508 : : //===========================================================================
509 : 198 : void FacetDataUtil::get_edges(
510 : : DLIList<CubitFacet *> &facet_list,
511 : : DLIList<CubitFacetEdge *> &edge_list )
512 : : {
513 : : int i, j;
514 : : CubitPoint *p0, *p1;
515 : : CubitFacet *facet_ptr;
516 : : CubitFacetEdge *edge_ptr;
517 [ + - ]: 198 : DLIList<CubitFacet *> adj_facet_list;
518 : :
519 : : // mark the edges and create any that are missing
520 [ + - ]: 198 : facet_list.reset(); //have to reset this list to that the CubitFacetEdgeData's
521 : : //get constructed correctly
522 [ + - ][ + + ]: 2574 : for ( i = 0; i < facet_list.size(); i++)
523 : : {
524 [ + - ]: 2376 : facet_ptr = facet_list.get_and_step();
525 [ + + ]: 9504 : for (j=0; j<3; j++) {
526 [ + - ]: 7128 : edge_ptr = facet_ptr->edge(j);
527 [ + + ]: 7128 : if (!(edge_ptr))
528 : : {
529 [ + - ]: 1694 : facet_ptr->get_edge_pts(j, p0, p1);
530 [ + - ][ + - ]: 1694 : edge_ptr = (CubitFacetEdge *) new CubitFacetEdgeData( p0, p1 );
531 : : }
532 [ + - ]: 7128 : edge_ptr->set_flag( 0 );
533 : : }
534 : : }
535 : :
536 : : // create a unique list of edges
537 : :
538 [ + - ][ + + ]: 2574 : for ( i = 0; i < facet_list.size(); i++)
539 : : {
540 [ + - ]: 2376 : facet_ptr = facet_list.get_and_step();
541 [ + + ]: 9504 : for (j=0; j<3; j++)
542 : : {
543 [ + - ]: 7128 : edge_ptr = facet_ptr->edge(j);
544 [ + - ][ + - ]: 7128 : if(edge_ptr && 0 == edge_ptr->get_flag())
[ + + ][ + + ]
545 : : {
546 [ + - ]: 3784 : edge_ptr->set_flag( 1 );
547 [ + - ]: 3784 : edge_list.append( edge_ptr );
548 : : }
549 : : }
550 : : }
551 : :
552 : : // reset the flags on the edges
553 : :
554 [ + - ][ + + ]: 2574 : for ( i = 0; i < facet_list.size(); i++)
555 : : {
556 [ + - ]: 2376 : facet_ptr = facet_list.get_and_step();
557 [ + + ]: 9504 : for (j=0; j<3; j++)
558 : : {
559 [ + - ]: 7128 : edge_ptr = facet_ptr->edge(j);
560 [ + - ]: 7128 : if( edge_ptr )
561 [ + - ]: 7128 : edge_ptr->set_flag( 0 );
562 : : }
563 [ + - ]: 198 : }
564 : 198 : }
565 : :
566 : : //============================================================================
567 : : // Function: quality
568 : : // Author: sjowen
569 : : // Description: this is the S.H. Lo metric for triangles that also takes into
570 : : // account a surface normal.
571 : : // Date: 2/2003
572 : : //============================================================================
573 : 0 : double FacetDataUtil::quality(CubitVector &c1, CubitVector &c2, CubitVector &c3,
574 : : CubitVector &surf_normal)
575 : : {
576 : : #define TWO_ROOT_THREE 3.46410161514
577 : : double area2, alpha;
578 : : double length1, length2, length3;
579 [ # # ][ # # ]: 0 : CubitVector edge1, edge2, edge3;
[ # # ]
580 : :
581 : : // create edges from the vertices
582 : :
583 [ # # ][ # # ]: 0 : edge1 = c3 - c2;
584 [ # # ][ # # ]: 0 : edge2 = c3 - c1;
585 [ # # ][ # # ]: 0 : edge3 = c2 - c1;
586 : :
587 : : // compute twice the area
588 : :
589 [ # # ]: 0 : CubitVector normal = edge3 * edge2;
590 [ # # ]: 0 : area2 = normal.length();
591 : :
592 [ # # ]: 0 : length1 = edge1.length_squared();
593 [ # # ]: 0 : length2 = edge2.length_squared();
594 [ # # ]: 0 : length3 = edge3.length_squared();
595 : :
596 : 0 : alpha = TWO_ROOT_THREE * area2 / (length1 + length2 + length3);
597 : :
598 : : // modify the alpha metric by the dot product of the triangle normal
599 : : // with the surface normal.
600 : :
601 [ # # ]: 0 : if (fabs(area2) < CUBIT_RESABS)
602 : 0 : alpha = 0.0;
603 : : else
604 : : {
605 [ # # ]: 0 : normal /= area2;
606 [ # # ]: 0 : double dot = normal % surf_normal;
607 : 0 : double penalty = pow(dot, 5);
608 : 0 : alpha *= penalty;
609 : : }
610 : :
611 : 0 : return alpha;
612 : : }
613 : :
614 : : //================================================================================
615 : : // Description: compares length of two edges.
616 : : // Notes: used for DLIList sort
617 : : // Author : Steve Owen
618 : : // Date : 9/27/2003
619 : : //================================================================================
620 : 0 : int FacetDataUtil::edge_compare(CubitFacetEdge *&ea, CubitFacetEdge *&eb)
621 : : {
622 : 0 : double la = ea->length();
623 : 0 : double lb = eb->length();
624 [ # # ]: 0 : if (la < lb)
625 : 0 : return -1;
626 [ # # ]: 0 : else if (lb < la)
627 : 0 : return 1;
628 : 0 : return 0;
629 : : }
630 : :
631 : : //================================================================================
632 : : // Description: collapse edges that are smaller than 10% the average length
633 : : // Notes: non_manifold_only=true indicates that only edges that are attached
634 : : // to more than two edges will be candidates for collapse.
635 : : // Author : Steve Owen
636 : : // Date : 9/27/2003
637 : : //================================================================================
638 : 11 : CubitStatus FacetDataUtil::collapse_short_edges(
639 : : DLIList<CubitFacet*> &facet_list,
640 : : CubitBoolean non_manifold_only)
641 : : {
642 : : #define COLLAPSE_TOLERANCE 0.3
643 : :
644 : 11 : CubitStatus rv = CUBIT_SUCCESS;
645 : 11 : int ncollapse = 0;
646 : :
647 : : // get the edges
648 : :
649 [ + - ]: 11 : DLIList <CubitFacetEdge *>edge_list;
650 [ + - ][ + - ]: 22 : DLIList <CubitPoint *> point_list;
651 [ + - ][ + - ]: 22 : DLIList <CubitFacet *> one_facet;
652 [ + - ][ + - ]: 22 : DLIList <CubitFacetEdge *> facet_edges;
653 [ + - ]: 11 : FacetDataUtil::get_edges( facet_list, edge_list );
654 [ + - ]: 11 : FacetDataUtil::get_points( facet_list, point_list );
655 : :
656 : : // determine average length and get the collapse threshold
657 : :
658 : : int ii, jj, kk;
659 : : double len;
660 : 11 : double tot_len = 0.0;
661 [ + - ][ + + ]: 99 : for(ii=0; ii<point_list.size(); ii++)
662 [ + - ][ + - ]: 88 : point_list.get_and_step()->marked(0);
663 : : CubitFacetEdge *edge;
664 [ + - ][ + + ]: 209 : for(ii=0; ii<edge_list.size(); ii++)
665 : : {
666 [ + - ]: 198 : edge = edge_list.get_and_step();
667 [ + - ]: 198 : len = edge->length();
668 : 198 : tot_len += len;
669 [ + - ]: 198 : edge->marked(0);
670 : : }
671 [ + - ]: 11 : len = tot_len / edge_list.size();
672 : 11 : double collapse_tol = COLLAPSE_TOLERANCE * len;
673 : : //check(edge_list, facet_list);
674 : :
675 : : // get a list of the edges that qualify for collapse
676 : :
677 [ + - ][ + - ]: 22 : DLIList<CubitFacetEdge *>short_edges;
678 [ + - ][ + + ]: 209 : for(ii=0; ii<edge_list.size(); ii++)
679 : : {
680 [ + - ]: 198 : edge = edge_list.get_and_step();
681 [ + - ][ + - ]: 198 : if (non_manifold_only && edge->num_adj_facets() == 2)
[ + - ][ + - ]
682 : 198 : continue;
683 [ # # ]: 0 : len = edge->length();
684 [ # # ]: 0 : if (len < collapse_tol)
685 : : {
686 [ # # ]: 0 : short_edges.append(edge);
687 : : }
688 : : }
689 : : //sort them
690 : :
691 [ + - ]: 11 : short_edges.sort(FacetDataUtil::edge_compare);
692 : :
693 : : // main loop
694 : :
695 [ + - ]: 11 : int nedges = short_edges.size();
696 [ - + ]: 11 : for(ii=0; ii<nedges; ii++)
697 : : {
698 [ # # ]: 0 : edge = short_edges.get_and_step();
699 [ # # ][ # # ]: 0 : if (edge->marked() == 1)
700 : 0 : continue;
701 [ # # ]: 0 : len = edge->length();
702 : :
703 : 0 : bool collapse = true;
704 : : CubitPoint *pt[2];
705 : : CubitPoint *pt1, *pt2;
706 [ # # ]: 0 : pt[0] = edge->point(0);
707 [ # # ]: 0 : pt[1] = edge->point(1);
708 : :
709 : : // compute a new candidate location for the merged points
710 : :
711 [ # # ]: 0 : CubitVector new_location;
712 [ # # ][ # # ]: 0 : CubitVector cpt[2];
713 [ # # ][ # # ]: 0 : cpt[0] = pt[0]->coordinates();
714 [ # # ][ # # ]: 0 : cpt[1] = pt[1]->coordinates();
715 [ # # ][ # # ]: 0 : new_location = (cpt[0] + cpt[1]) * 0.5;
[ # # ]
716 : :
717 [ # # ][ # # ]: 0 : for (jj=0; jj<2 && collapse; jj++)
718 : : {
719 [ # # ]: 0 : DLIList<CubitFacet *> adjfacets;
720 [ # # ]: 0 : pt[jj]->facets( adjfacets );
721 : :
722 : : // check all facets adjacent this point that don't contain the edge
723 : : // to make sure the resulting facet will be valid (we aren't inverting anything)
724 : :
725 [ # # ][ # # ]: 0 : for(kk=0; kk<adjfacets.size() && collapse; kk++)
[ # # ][ # # ]
726 : : {
727 [ # # ]: 0 : CubitFacet *facet = adjfacets.get_and_step();
728 [ # # ]: 0 : pt1 = facet->next_node( pt[jj] );
729 [ # # ]: 0 : pt2 = facet->next_node( pt1 );
730 [ # # ]: 0 : int eidx = facet->edge_index( edge );
731 [ # # ]: 0 : if (eidx < 0)
732 : : // don't check facets that have the current edge - they'll be deleted anyway
733 : : {
734 : :
735 [ # # ]: 0 : CubitVector cpt1 = pt1->coordinates();
736 [ # # ]: 0 : CubitVector cpt2 = pt2->coordinates();
737 [ # # ]: 0 : CubitVector norm = facet->normal();
738 [ # # ]: 0 : double q0 = FacetDataUtil::quality(cpt[jj], cpt1, cpt2, norm);
739 [ # # ]: 0 : double q1 = FacetDataUtil::quality(new_location, cpt1, cpt2, norm);
740 [ # # ][ # # ]: 0 : if (!(q1 > 0.0 || q1 > q0))
741 : : {
742 : 0 : collapse = false;
743 : : }
744 : : }
745 : : }
746 [ # # ]: 0 : }
747 : :
748 [ # # ]: 0 : if (collapse)
749 : : {
750 [ # # ]: 0 : DLIList<CubitFacet *> new_facets;
751 : 0 : CubitPoint *collpt = pt[0];
752 : 0 : CubitPoint *delpt = pt[1];
753 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *> adjedges;
754 : : CubitFacetEdge *adjedge;
755 [ # # ]: 0 : delpt->edges( adjedges );
756 : : CubitFacet *facet;
757 : :
758 : :
759 : 0 : int mydebug = 0;
760 [ # # ]: 0 : if (mydebug)
761 : : {
762 [ # # ]: 0 : dcolor(4);
763 [ # # ]: 0 : draw_edge( edge );
764 [ # # ]: 0 : dview();
765 : : }
766 : :
767 [ # # ]: 0 : collpt->set(new_location);
768 : :
769 : : // delete all facets adjacent to the delpt.
770 : :
771 [ # # ][ # # ]: 0 : DLIList<CubitFacet *> adjfacets;
772 [ # # ]: 0 : delpt->facets( adjfacets );
773 : :
774 [ # # ][ # # ]: 0 : for(jj=0; jj<adjfacets.size(); jj++)
775 : : {
776 [ # # ]: 0 : facet = adjfacets.get_and_step();
777 [ # # ]: 0 : pt1 = facet->next_node(delpt);
778 [ # # ]: 0 : pt2 = facet->next_node(pt1);
779 [ # # ]: 0 : int eidx = facet->edge_index( edge );
780 [ # # ]: 0 : facet_list.move_to(facet);
781 [ # # ][ # # ]: 0 : delete facet;
782 : 0 : facet = NULL;
783 [ # # ]: 0 : if (eidx >= 0)
784 : : {
785 : : //if this facet is adjecnt the edge, then just remove it from the list
786 [ # # ]: 0 : facet_list.extract();
787 : : }
788 : : else
789 : : {
790 : :
791 : : // get or create edges as needed
792 : :
793 : : CubitFacetEdge *e[3];
794 [ # # ]: 0 : e[0] = collpt->get_edge( pt1 );
795 [ # # ]: 0 : e[1] = pt1->get_edge( pt2 );
796 [ # # ]: 0 : e[2] = pt2->get_edge( collpt );
797 [ # # ]: 0 : for (kk=0; kk<3; kk++)
798 : : {
799 [ # # ]: 0 : if (!(e[kk]))
800 : : {
801 [ # # # # ]: 0 : switch (kk)
802 : : {
803 [ # # ][ # # ]: 0 : case 0: e[kk] = (CubitFacetEdge *) new CubitFacetEdgeData( collpt, pt1 ); break;
804 [ # # ][ # # ]: 0 : case 1: e[kk] = (CubitFacetEdge *) new CubitFacetEdgeData( pt1, pt2 ); break;
805 [ # # ][ # # ]: 0 : case 2: e[kk] = (CubitFacetEdge *) new CubitFacetEdgeData( pt2, collpt ); break;
806 : : }
807 [ # # ]: 0 : edge_list.append( e[kk] );
808 : : }
809 : : }
810 : :
811 : : // create a new facet with the points from the old facet and the collpt
812 : :
813 [ # # ][ # # ]: 0 : facet = new CubitFacetData( e[0], e[1], e[2] );
814 [ # # ]: 0 : new_facets.append(facet);
815 : :
816 : : // create edges on the facet
817 : :
818 [ # # ]: 0 : facet_list.change_to( facet );
819 : : }
820 : : }
821 : :
822 [ # # ][ # # ]: 0 : for(jj=0; jj<adjedges.size(); jj++)
823 : : {
824 [ # # ]: 0 : adjedge = adjedges.get_and_step();
825 [ # # ]: 0 : adjedge->marked(1); // mark it for deletion later
826 [ # # ][ # # ]: 0 : assert(adjedge->num_adj_facets() == 0);
827 : : }
828 [ # # ][ # # ]: 0 : assert(delpt->num_adj_facets() == 0);
829 [ # # ]: 0 : ncollapse++;
830 : :
831 : : //check(edge_list, facet_list);
832 : : }
833 : :
834 : : }
835 : :
836 [ + - ][ + - ]: 11 : PRINT_INFO("Collapsed %d short edges in triangulation\n", ncollapse);
[ + - ][ + - ]
837 : :
838 : : // delete points and edges that aren't used
839 [ - + ]: 11 : if (ncollapse > 0)
840 : : {
841 [ # # ][ # # ]: 0 : for(ii=0; ii<edge_list.size(); ii++)
842 : : {
843 [ # # ]: 0 : edge = edge_list.get_and_step();
844 [ # # ][ # # ]: 0 : if (edge->marked())
845 : : {
846 [ # # ][ # # ]: 0 : assert(edge->num_adj_facets() == 0);
847 [ # # ][ # # ]: 0 : delete edge;
848 : : }
849 [ # # ][ # # ]: 0 : else if(edge->num_adj_facets() < 2){
850 [ # # ][ # # ]: 0 : PRINT_ERROR("Unexpected result while collapsing an edge.\n");
[ # # ][ # # ]
851 : 0 : return CUBIT_FAILURE;
852 : : //assert(edge->num_adj_facets() >= 2);
853 : : }
854 : : }
855 : : CubitPoint *point;
856 [ # # ][ # # ]: 0 : for(ii=0; ii<point_list.size(); ii++)
857 : : {
858 [ # # ]: 0 : point = point_list.get_and_step();
859 [ # # ][ # # ]: 0 : if (point->num_adj_facets() == 0)
860 : : {
861 [ # # ][ # # ]: 0 : delete point;
862 : : }
863 : : }
864 : : }
865 : :
866 [ + - ]: 22 : return rv;
867 : : }
868 : :
869 : : //===========================================================================
870 : : // Function: check
871 : : // Purpose: debugging only
872 : : // Date: 10/2003
873 : : // Author: sjowen
874 : : //===========================================================================
875 : 0 : void FacetDataUtil::check(DLIList<CubitFacetEdge *> &edge_list,
876 : : DLIList<CubitFacet *> &facet_list)
877 : : {
878 : :
879 [ # # ]: 0 : CubitBox box;
880 : : CubitFacetEdge *edge;
881 : : int ii;
882 : 0 : int nedges = 0;
883 [ # # ][ # # ]: 0 : for(ii=0; ii<edge_list.size(); ii++)
884 : : {
885 [ # # ]: 0 : edge = edge_list.get_and_step();
886 [ # # ][ # # ]: 0 : if (edge->marked() == 1)
887 : 0 : continue;
888 [ # # ]: 0 : int nadj = edge->num_adj_facets();
889 [ # # ]: 0 : if (nadj <= 1)
890 : : {
891 [ # # ][ # # ]: 0 : CubitVector v0 = edge->point(0)->coordinates();
892 [ # # ][ # # ]: 0 : CubitVector v1 = edge->point(1)->coordinates();
893 [ # # ][ # # ]: 0 : CubitVector min(CUBIT_MIN(v0.x(), v1.x()), CUBIT_MIN(v0.y(), v1.y()), CUBIT_MIN(v0.z(), v1.z()));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
894 [ # # ][ # # ]: 0 : CubitVector max(CUBIT_MAX(v0.x(), v1.x()), CUBIT_MAX(v0.y(), v1.y()), CUBIT_MAX(v0.z(), v1.z()));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
895 [ # # ]: 0 : CubitBox ebox(min, max);
896 [ # # ]: 0 : if (nedges == 0)
897 [ # # ]: 0 : box.reset(min, max);
898 : : else
899 [ # # ]: 0 : box |= ebox;
900 : : // dcolor(3);
901 : : // dfldraw(facet_list);
902 [ # # ]: 0 : dcolor(4);
903 [ # # ]: 0 : dedraw(edge);
904 [ # # ]: 0 : dpoint(v0);
905 [ # # ]: 0 : dpoint(v1);
906 [ # # ]: 0 : nedges++;
907 : : }
908 : : }
909 [ # # ]: 0 : if (nedges > 0)
910 : : {
911 [ # # ]: 0 : dzoom(box);
912 [ # # ]: 0 : dview();
913 [ # # ]: 0 : }
914 : 0 : }
915 : :
916 : : //===========================================================================
917 : : // Function: draw_edge
918 : : // Purpose:
919 : : // Date: 10/2003
920 : : // Author: sjowen
921 : : //===========================================================================
922 : 0 : void FacetDataUtil::draw_edge( CubitFacetEdge *edge )
923 : : {
924 [ # # ]: 0 : DLIList<CubitFacet *>adjfacets;
925 [ # # ]: 0 : edge->facets( adjfacets );
926 : : int ii,jj;
927 : : CubitFacet *facet;
928 [ # # ][ # # ]: 0 : CubitBox box;
929 [ # # ]: 0 : for(jj=0; jj<2; jj++)
930 : : {
931 [ # # ]: 0 : CubitPoint *p = edge->point(jj);
932 [ # # ]: 0 : adjfacets.clean_out();
933 [ # # ]: 0 : p->facets( adjfacets );
934 [ # # ][ # # ]: 0 : for(ii=0; ii<adjfacets.size(); ii++)
935 : : {
936 [ # # ]: 0 : facet = adjfacets.get_and_step();
937 [ # # ][ # # ]: 0 : if (jj==0 && ii==0)
938 [ # # ][ # # ]: 0 : box = facet->bounding_box();
939 : : else
940 [ # # ][ # # ]: 0 : box |= facet->bounding_box();
941 [ # # ]: 0 : dfdraw(facet);
942 : : }
943 [ # # ][ # # ]: 0 : dpoint(p->coordinates());
944 [ # # ][ # # ]: 0 : dpoint(p->coordinates());
945 : : }
946 [ # # ][ # # ]: 0 : dzoom(box);
947 : 0 : }
948 : :
949 : : //=============================================================================
950 : : //Function: is_point_in_polyhedron
951 : : //Description: point-in-polyhedron test; polyhedron must be water-tight,
952 : : // manifold, triangle-tiled. Casts a random ray in positive
953 : : // x, y, and z. Counts crossings. Starts over with a recast
954 : : // if a triangle is perpendicular to the ray or if the ray
955 : : // hits a triangle edge or vertex.
956 : : // Also tests for point on polyhedron.
957 : : //Author: John Fowler
958 : : //Date: 9/30/03
959 : : //=============================================================================
960 : 0 : CubitStatus FacetDataUtil::is_point_in_polyhedron(
961 : : DLIList<CubitFacet *> &tfacet_list,
962 : : CubitVector &point_coords,
963 : : CubitPointContainment &is_point_in)
964 : : {
965 : : unsigned int i;
966 [ # # ]: 0 : CubitBox bbox;
967 : : CubitFacet *facet;
968 [ # # ][ # # ]: 0 : CubitVector bbox_min, bbox_max, ray;
[ # # ]
969 : : CubitStatus status;
970 : : CubitPointContainment pt_status;
971 : : bool is_outside, recast, is_on_plane;
972 : : double xpt, ypt, zpt;
973 : : double rayx, rayy, rayz;
974 : : int number_of_recasts;
975 : :
976 [ # # ]: 0 : point_coords.get_xyz(xpt, ypt, zpt);
977 : 0 : recast = true;
978 : 0 : number_of_recasts = 0;
979 [ # # ][ # # ]: 0 : while ( (recast == true) && (number_of_recasts < 10) ) {
980 : 0 : recast = false;
981 : 0 : is_outside = true;
982 [ # # ]: 0 : random_positive_ray(ray); // a random positively-pointing ray
983 [ # # ]: 0 : ray.get_xyz(rayx,rayy,rayz);
984 [ # # ][ # # ]: 0 : for ( i = tfacet_list.size(); i > 0; i-- ) {
985 [ # # ]: 0 : facet = tfacet_list.get_and_step();
986 [ # # ][ # # ]: 0 : bbox = facet->bounding_box();
987 [ # # ][ # # ]: 0 : bbox_min = bbox.minimum();
988 [ # # ][ # # ]: 0 : bbox_max = bbox.maximum();
989 : : // Because the ray is positive in direction, discard bounding boxes
990 : : // that are entirely < the starting point.
991 [ # # ][ # # ]: 0 : if ( (xpt-bbox_max.x()) > CUBIT_RESABS ) continue;
992 [ # # ][ # # ]: 0 : if ( (ypt-bbox_max.y()) > CUBIT_RESABS ) continue;
993 [ # # ][ # # ]: 0 : if ( (zpt-bbox_max.z()) > CUBIT_RESABS ) continue;
994 [ # # ][ # # ]: 0 : if ( ray_intersects_boundingbox(point_coords,ray,bbox) == false )
995 : 0 : continue;
996 [ # # ][ # # ]: 0 : CubitPlane plane = facet->plane();
997 [ # # ][ # # ]: 0 : CubitVector normal = plane.normal();
998 : : double xnorm, ynorm, znorm;
999 [ # # ]: 0 : xnorm = normal.x();
1000 [ # # ]: 0 : ynorm = normal.y();
1001 [ # # ]: 0 : znorm = normal.z();
1002 : : // Intersect the ray with the facet plane. If ray is perpendicular to
1003 : : // facet plane and point is on it, recast the ray and try again.
1004 : 0 : double denominator = rayx*xnorm + rayy*ynorm + rayz*znorm;
1005 [ # # ]: 0 : double distanc = xnorm*xpt + ynorm*ypt + znorm*zpt + plane.coefficient();
1006 [ # # ]: 0 : if ( fabs(denominator) < GEOMETRY_RESABS )
1007 : : {
1008 [ # # ]: 0 : if ( fabs(distanc) < GEOMETRY_RESABS ) {
1009 : 0 : recast = true;
1010 : 0 : break;
1011 : 0 : } else continue; // point is not on plane and ray is parallel to plane
1012 : : }
1013 : : double t, xintpt, yintpt, zintpt;
1014 : 0 : t = -distanc;
1015 : 0 : t /= denominator;
1016 [ # # ]: 0 : if ( fabs(t) < GEOMETRY_RESABS ) is_on_plane = true;
1017 [ # # ]: 0 : else if ( t < GEOMETRY_RESABS ) continue;
1018 : 0 : else is_on_plane = false;
1019 : 0 : xintpt = xpt + t*rayx;
1020 : 0 : yintpt = ypt + t*rayy;
1021 : 0 : zintpt = zpt + t*rayz;
1022 : : CubitPoint *pt;
1023 : : // We need to project the triangle onto 2D. Use the smaller components of
1024 : : // the normal to detemine a good projection.
1025 : : double x0, y0, x1, y1, x2, y2;
1026 [ # # ][ # # ]: 0 : if ( (fabs(xnorm) >= fabs(ynorm)) && (fabs(xnorm) >= fabs(znorm)) ){ // Use y,z
1027 [ # # ]: 0 : pt = facet->point(0);
1028 [ # # ][ # # ]: 0 : x0 = pt->y(); y0 = pt->z();
1029 [ # # ]: 0 : pt = facet->point(1);
1030 [ # # ][ # # ]: 0 : x1 = pt->y(); y1 = pt->z();
1031 [ # # ]: 0 : pt = facet->point(2);
1032 [ # # ][ # # ]: 0 : x2 = pt->y(); y2 = pt->z();
1033 [ # # ]: 0 : status = pt_in_tri_2d(yintpt,zintpt,x0,y0,x1,y1,x2,y2,pt_status);
1034 [ # # ]: 0 : } else if (fabs(ynorm) >= fabs(znorm)) { // Use z,x
1035 [ # # ]: 0 : pt = facet->point(0);
1036 [ # # ][ # # ]: 0 : x0 = pt->x(); y0 = pt->z();
1037 [ # # ]: 0 : pt = facet->point(1);
1038 [ # # ][ # # ]: 0 : x1 = pt->x(); y1 = pt->z();
1039 [ # # ]: 0 : pt = facet->point(2);
1040 [ # # ][ # # ]: 0 : x2 = pt->x(); y2 = pt->z();
1041 [ # # ]: 0 : status = pt_in_tri_2d(xintpt,zintpt,x0,y0,x1,y1,x2,y2,pt_status);
1042 : : } else { // Use x,y
1043 [ # # ]: 0 : pt = facet->point(0);
1044 [ # # ][ # # ]: 0 : x0 = pt->x(); y0 = pt->y();
1045 [ # # ]: 0 : pt = facet->point(1);
1046 [ # # ][ # # ]: 0 : x1 = pt->x(); y1 = pt->y();
1047 [ # # ]: 0 : pt = facet->point(2);
1048 [ # # ][ # # ]: 0 : x2 = pt->x(); y2 = pt->y();
1049 [ # # ]: 0 : status = pt_in_tri_2d(xintpt,yintpt,x0,y0,x1,y1,x2,y2,pt_status);
1050 : : }
1051 : :
1052 [ # # ]: 0 : if ( status == CUBIT_FAILURE ) {
1053 : 0 : recast = true;
1054 : 0 : break;
1055 : : }
1056 [ # # ]: 0 : if ( pt_status == CUBIT_PNT_OUTSIDE ) continue;
1057 : : // Is the point on the triangle?
1058 [ # # ]: 0 : if ( is_on_plane == true ) {
1059 : 0 : is_point_in = CUBIT_PNT_BOUNDARY;
1060 : 0 : return CUBIT_SUCCESS;
1061 : : }
1062 [ # # ]: 0 : if ( pt_status == CUBIT_PNT_INSIDE ) is_outside = ! is_outside;
1063 [ # # ]: 0 : else if ( pt_status == CUBIT_PNT_BOUNDARY ) {
1064 : 0 : recast = true;
1065 : 0 : break;
1066 : : }
1067 : : }
1068 [ # # ]: 0 : if ( recast == true ) {
1069 [ # # ]: 0 : tfacet_list.reset();
1070 : 0 : number_of_recasts += 1;
1071 : : }
1072 : : }
1073 [ # # ]: 0 : if ( recast == true ) {
1074 [ # # ][ # # ]: 0 : PRINT_ERROR("Number of recasts in point-in-polygon exceeded 10.\n");
[ # # ][ # # ]
1075 : 0 : return CUBIT_FAILURE;
1076 : : }
1077 [ # # ]: 0 : if ( is_outside == false ) is_point_in = CUBIT_PNT_INSIDE;
1078 : 0 : else is_point_in = CUBIT_PNT_OUTSIDE;
1079 : :
1080 [ # # ]: 0 : return CUBIT_SUCCESS;
1081 : : }
1082 : :
1083 : : //===========================================================================
1084 : : // Function: pt_in_tri_2d
1085 : : // Purpose:
1086 : : // Date: 10/2003
1087 : : // Author: John Fowler
1088 : : //===========================================================================
1089 : 0 : CubitStatus FacetDataUtil::pt_in_tri_2d(double xpt, double ypt,
1090 : : double x0, double y0,
1091 : : double x1, double y1,
1092 : : double x2, double y2,
1093 : : CubitPointContainment &is_point_in)
1094 : : {
1095 : : // From Schneider & Eberly, "Geometric Tools for COmputer Graphics",
1096 : : // Chap. 13.3.1. If triangle is needle-thin, CUBIT_FAILURE might be
1097 : : // returned, in wich case is_point_in is undefined.
1098 : :
1099 : : double c0, c1, c2;
1100 : : double e0x, e1x, e2x, e0y, e1y, e2y;
1101 : : double n0x, n1x, n2x, n0y, n1y, n2y;
1102 : : double denom0, denom1, denom2;
1103 : :
1104 : 0 : e0x = x1 - x0; e0y = y1 - y0;
1105 : 0 : e1x = x2 - x1; e1y = y2 - y1;
1106 : 0 : e2x = x0 - x2; e2y = y0 - y2;
1107 : 0 : n0x = e0y; n0y = -e0x;
1108 : 0 : n1x = e1y; n1y = -e1x;
1109 : 0 : n2x = e2y; n2y = -e2x;
1110 : 0 : denom0 = n1x*e0x + n1y*e0y;
1111 [ # # ]: 0 : if ( fabs(denom0) < CUBIT_RESABS ) {
1112 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
1113 : 0 : return CUBIT_FAILURE;
1114 : : }
1115 : 0 : denom1 = n2x*e1x + n2y*e1y;
1116 [ # # ]: 0 : if ( fabs(denom1) < CUBIT_RESABS ) {
1117 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
1118 : 0 : return CUBIT_FAILURE;
1119 : : }
1120 : 0 : denom2 = n0x*e2x + n0y*e2y;
1121 [ # # ]: 0 : if ( fabs(denom2) < CUBIT_RESABS ) {
1122 [ # # ][ # # ]: 0 : PRINT_ERROR("Failure in pt_in_tri_2d; needle-thin triangle encountered.\n");
1123 : 0 : return CUBIT_FAILURE;
1124 : : }
1125 : :
1126 : 0 : c0 = -( n1x*(xpt-x1) + n1y*(ypt-y1) )/denom0;
1127 : 0 : c1 = -( n2x*(xpt-x2) + n2y*(ypt-y2) )/denom1;
1128 : 0 : c2 = -( n0x*(xpt-x0) + n0y*(ypt-y0) )/denom2;
1129 : :
1130 [ # # ][ # # ]: 0 : if ( (c0 > 0.0) && (c1 > 0.0) && (c2 > 0.0) ) is_point_in = CUBIT_PNT_INSIDE;
[ # # ]
1131 [ # # ][ # # ]: 0 : else if ( (c0 < 0.0) || (c1 < 0.0) || (c2 < 0.0) ) is_point_in = CUBIT_PNT_OUTSIDE;
[ # # ]
1132 : 0 : else is_point_in = CUBIT_PNT_BOUNDARY;
1133 : :
1134 : 0 : return CUBIT_SUCCESS;
1135 : :
1136 : : }
1137 : :
1138 : : //===========================================================================
1139 : : // Function: random_positive_ray
1140 : : // Purpose:
1141 : : // Date: 10/2003
1142 : : // Author: John Fowler
1143 : : //===========================================================================
1144 : 0 : void FacetDataUtil::random_positive_ray(CubitVector& ray)
1145 : : {
1146 : : double temp;
1147 : 0 : double rayx = 0.0, rayy = 0.0, rayz = 0.0;
1148 : :
1149 : 0 : temp = 0.0;
1150 [ # # ]: 0 : while ( temp < 1.e-6 ) {
1151 : 0 : rayx = (double(rand())/(RAND_MAX+1.0));
1152 : 0 : rayy = (double(rand())/(RAND_MAX+1.0));
1153 : 0 : rayz = (double(rand())/(RAND_MAX+1.0));
1154 : 0 : temp = sqrt(rayx*rayx + rayy*rayy + rayz*rayz);
1155 : : }
1156 : 0 : rayx /= temp;
1157 : 0 : rayy /= temp;
1158 : 0 : rayz /= temp;
1159 : :
1160 : 0 : ray.set(rayx,rayy,rayz);
1161 : :
1162 : 0 : }
1163 : :
1164 : : //===========================================================================
1165 : : // Function: ray_intersects_boundingbox
1166 : : // Purpose:
1167 : : // Date: 10/2003
1168 : : // Author: John Fowler
1169 : : //===========================================================================
1170 : 0 : bool FacetDataUtil::ray_intersects_boundingbox(CubitVector& point, CubitVector& ray, const CubitBox& bbox)
1171 : : {
1172 : : double t, xtest, ytest, ztest;
1173 : : double xdir, ydir, zdir, xpt, ypt, zpt, xmin, ymin, zmin, xmax, ymax, zmax;
1174 [ # # ][ # # ]: 0 : CubitVector bbox_min, bbox_max;
1175 : :
1176 [ # # ]: 0 : point.get_xyz(xpt,ypt,zpt);
1177 [ # # ]: 0 : ray.get_xyz(xdir,ydir,zdir);
1178 [ # # ][ # # ]: 0 : bbox_min = bbox.minimum();
1179 [ # # ][ # # ]: 0 : bbox_max = bbox.maximum();
1180 [ # # ]: 0 : bbox_min.get_xyz(xmin,ymin,zmin);
1181 [ # # ]: 0 : bbox_max.get_xyz(xmax,ymax,zmax);
1182 : 0 : xmin -= GEOMETRY_RESABS; // So we don't miss any.
1183 : 0 : ymin -= GEOMETRY_RESABS;
1184 : 0 : zmin -= GEOMETRY_RESABS;
1185 : 0 : xmax += GEOMETRY_RESABS;
1186 : 0 : ymax += GEOMETRY_RESABS;
1187 : 0 : zmax += GEOMETRY_RESABS;
1188 : :
1189 : : // Notice that we are only interested in bboxes on the non-negative side of t.
1190 [ # # ]: 0 : if ( fabs(xdir) > CUBIT_RESABS ) {
1191 : : // test xmin plane
1192 : 0 : t = (xmin - xpt)/xdir;
1193 [ # # ]: 0 : if ( t >= 0.0 ) {
1194 : 0 : ytest = ypt + t*ydir;
1195 : 0 : ztest = zpt + t*zdir;
1196 [ # # ][ # # ]: 0 : if ( (ytest >= ymin) && (ytest <= ymax) && (ztest >= zmin) && (ztest <= zmax) )
[ # # ][ # # ]
1197 : 0 : return true;
1198 : : }
1199 : : // test xmax plane
1200 : 0 : t = (xmax - xpt)/xdir;
1201 [ # # ]: 0 : if ( t > 0.0 ) {
1202 : 0 : ytest = ypt + t*ydir;
1203 : 0 : ztest = zpt + t*zdir;
1204 [ # # ][ # # ]: 0 : if ( (ytest >= ymin) && (ytest <= ymax) && (ztest >= zmin) && (ztest <= zmax) )
[ # # ][ # # ]
1205 : 0 : return true;
1206 : : }
1207 : : }
1208 [ # # ]: 0 : if ( fabs(ydir) > CUBIT_RESABS ) {
1209 : : // test ymin plane
1210 : 0 : t = (ymin - ypt)/ydir;
1211 [ # # ]: 0 : if ( t >= 0.0 ) {
1212 : 0 : xtest = xpt + t*xdir;
1213 : 0 : ztest = zpt + t*zdir;
1214 [ # # ][ # # ]: 0 : if ( (xtest >= xmin) && (xtest <= xmax) && (ztest >= zmin) && (ztest <= zmax) )
[ # # ][ # # ]
1215 : 0 : return true;
1216 : : }
1217 : : // test ymax plane
1218 : :
1219 : 0 : t = (ymax - ypt)/ydir;
1220 [ # # ]: 0 : if ( t > 0.0 ) {
1221 : 0 : xtest = xpt + t*xdir;
1222 : 0 : ztest = zpt + t*zdir;
1223 [ # # ][ # # ]: 0 : if ( (xtest >= xmin) && (xtest <= xmax) && (ztest >= zmin) && (ztest <= zmax) )
[ # # ][ # # ]
1224 : 0 : return true;
1225 : : }
1226 : : }
1227 [ # # ]: 0 : if ( fabs(zdir) > CUBIT_RESABS ) {
1228 : : // test zmin plane
1229 : 0 : t = (zmin - zpt)/zdir;
1230 [ # # ]: 0 : if ( t > 0.0 ) {
1231 : 0 : xtest = xpt + t*xdir;
1232 : 0 : ytest = ypt + t*ydir;
1233 [ # # ][ # # ]: 0 : if ( (xtest >= xmin) && (xtest <= xmax) && (ytest >= ymin) && (ytest <= ymax) )
[ # # ][ # # ]
1234 : 0 : return true;
1235 : : }
1236 : : // test zmax plane
1237 : 0 : t = (zmax - zpt)/zdir;
1238 [ # # ]: 0 : if ( t > 0.0 ) {
1239 : 0 : xtest = xpt + t*xdir;
1240 : 0 : ytest = ypt + t*ydir;
1241 [ # # ][ # # ]: 0 : if ( (xtest >= xmin) && (xtest <= xmax) && (ytest >= ymin) && (ytest <= ymax) )
[ # # ][ # # ]
1242 : 0 : return true;
1243 : : }
1244 : : }
1245 : :
1246 : 0 : return false;
1247 : : }
1248 : :
1249 : : //===========================================================================
1250 : : // Function: write_facets
1251 : : // Purpose: write a list of facets to a cubit facet file
1252 : : // Date: 11/28/2002
1253 : : // Author: sjowen
1254 : : //===========================================================================
1255 : 0 : CubitStatus FacetDataUtil::write_facets( const char *file_name, DLIList<CubitFacet *> &facet_list)
1256 : : {
1257 [ # # ]: 0 : FILE *fp = fopen(file_name, "w");
1258 [ # # ]: 0 : if (!fp)
1259 : : {
1260 [ # # ][ # # ]: 0 : PRINT_ERROR("Couldn't open file %s for writing a facet file.\n", file_name);
[ # # ][ # # ]
1261 : 0 : return CUBIT_FAILURE;
1262 : : }
1263 : :
1264 [ # # ]: 0 : DLIList<CubitPoint*> point_list;
1265 [ # # ]: 0 : get_points(facet_list, point_list);
1266 : :
1267 [ # # ][ # # ]: 0 : fprintf(fp, "%d %d\n", point_list.size(), facet_list.size());
[ # # ]
1268 : : CubitPoint *pt;
1269 : : int ii;
1270 [ # # ][ # # ]: 0 : for (ii=1; ii<=point_list.size(); ii++)
1271 : : {
1272 [ # # ]: 0 : pt = point_list.get_and_step();
1273 [ # # ]: 0 : pt->set_id(ii);
1274 [ # # ][ # # ]: 0 : fprintf(fp, "%d %f %f %f\n", ii, pt->x(), pt->y(), pt->z() );
[ # # ][ # # ]
1275 : : }
1276 : :
1277 : : CubitFacet *facet;
1278 [ # # ][ # # ]: 0 : for (ii=1; ii<=facet_list.size(); ii++)
1279 : : {
1280 [ # # ]: 0 : facet = facet_list.get_and_step();
1281 [ # # ]: 0 : fprintf(fp, "%d %d %d %d\n", ii, facet->point(0)->id(),
1282 [ # # ][ # # ]: 0 : facet->point(1)->id(), facet->point(2)->id());
[ # # ][ # # ]
[ # # ][ # # ]
1283 : : }
1284 : :
1285 [ # # ]: 0 : fclose(fp);
1286 [ # # ]: 0 : return CUBIT_SUCCESS;
1287 : :
1288 : : }
1289 : :
1290 : : //=============================================================================
1291 : : //Function: split_into_shells (PUBLIC)
1292 : : //Description: split this set of facets into multiple surfaces where there are
1293 : : // may be disjoint regions.
1294 : : //Author: sjowen
1295 : : //Date: 8/27/2003
1296 : : //=============================================================================
1297 : 11 : CubitStatus FacetDataUtil::split_into_shells(
1298 : : DLIList<CubitFacet *> &tfacet_list,
1299 : : DLIList<CubitQuadFacet *> &qfacet_list,
1300 : : DLIList<DLIList<CubitFacet *> *> &shell_list,
1301 : : CubitBoolean &is_water_tight)
1302 : : {
1303 : 11 : CubitStatus stat = CUBIT_SUCCESS;
1304 : : //need to init this variable otherwise the caller must have.
1305 : 11 : is_water_tight = CUBIT_TRUE;
1306 : :
1307 : : // combine the quads and tri lists
1308 : :
1309 [ + - ]: 11 : DLIList <CubitFacet *> facet_list = tfacet_list;
1310 : : int ii;
1311 : : CubitQuadFacet *qfacet_ptr;
1312 [ + - ][ - + ]: 11 : for (ii=0; ii<qfacet_list.size(); ii++)
1313 : : {
1314 [ # # ]: 0 : qfacet_ptr = qfacet_list.get_and_step();
1315 [ # # ][ # # ]: 0 : facet_list.append(qfacet_ptr->get_tri_facet( 0 ));
1316 [ # # ][ # # ]: 0 : facet_list.append(qfacet_ptr->get_tri_facet( 1 ));
1317 : : }
1318 : :
1319 : : // mark all the facets to begin with
1320 : :
1321 : : CubitFacet *facet_ptr;
1322 [ + - ][ + + ]: 143 : for (ii=0; ii<facet_list.size(); ii++)
1323 : : {
1324 [ + - ]: 132 : facet_ptr = facet_list.get_and_step();
1325 [ + - ]: 132 : facet_ptr->marked( 1 );
1326 : : }
1327 : :
1328 : : // populate the facet edge lists
1329 : :
1330 [ + - ][ + - ]: 22 : DLIList<CubitFacetEdge *> edge_list;
1331 [ + - ]: 11 : FacetDataUtil::get_edges( facet_list, edge_list );
1332 : :
1333 : : // some debug stuff to draw the facets
1334 : :
1335 : 11 : int mydebug=0;
1336 [ - + ]: 11 : if (mydebug)
1337 : : {
1338 [ # # ][ # # ]: 0 : for (ii=0; ii<facet_list.size(); ii++)
1339 : : {
1340 [ # # ]: 0 : facet_ptr = facet_list.get_and_step();
1341 [ # # ]: 0 : GfxDebug::draw_facet(facet_ptr,CUBIT_YELLOW_INDEX);
1342 : : }
1343 [ # # ]: 0 : GfxDebug::flush();
1344 [ # # ]: 0 : GfxDebug::mouse_xforms();
1345 : : }
1346 : :
1347 : : // Go through the surfaceElemList and pull facets off one by one as we
1348 : : // determine which surface it belongs to. Continue until we have depleted
1349 : : // the list
1350 : :
1351 : : int jj;
1352 : 11 : int num_surfs_created = 0;
1353 [ + - ]: 11 : int num_facets_remaining = facet_list.size();
1354 [ + + ]: 22 : while( num_facets_remaining > 0)
1355 : : {
1356 : : // create a new shell to hold the face info
1357 : :
1358 : 11 : num_surfs_created++;
1359 [ + - ][ + - ]: 11 : DLIList<CubitFacet *> *shell_ptr = new DLIList<CubitFacet *>;
1360 : :
1361 : : // start with the first facet and create a list of all elements
1362 : : // attached to the facet
1363 : :
1364 : 11 : CubitBoolean shell_is_water_tight = CUBIT_TRUE;
1365 [ + - ]: 11 : CubitFacet *start_facet_ptr = facet_list.get_and_step();
1366 : : stat = get_adj_facets_on_shell( start_facet_ptr, shell_ptr,
1367 [ + - ]: 11 : shell_is_water_tight, mydebug );
1368 [ - + ]: 11 : if (stat != CUBIT_SUCCESS)
1369 : 0 : return stat;
1370 [ - + ]: 11 : if (!shell_is_water_tight)
1371 : 0 : is_water_tight = CUBIT_FALSE;
1372 : :
1373 [ + - ]: 11 : shell_list.append( shell_ptr );
1374 : :
1375 : : // remove the facets in this shell from the facet list
1376 : :
1377 [ + - ][ + + ]: 143 : for( jj = facet_list.size(); jj > 0; jj-- )
1378 : : {
1379 [ + - ]: 132 : facet_ptr = facet_list.get();
1380 [ + - ][ + - ]: 132 : if (facet_ptr->marked() == 0)
1381 : : {
1382 [ + - ]: 132 : facet_list.change_to( NULL );
1383 : : }
1384 [ + - ]: 132 : facet_list.step();
1385 : : }
1386 [ + - ]: 11 : facet_list.remove_all_with_value( NULL );
1387 [ + - ]: 11 : num_facets_remaining = facet_list.size();
1388 : : }
1389 : :
1390 [ + - ]: 22 : return CUBIT_SUCCESS;
1391 : : }
1392 : :
1393 : : //=============================================================================
1394 : : //Function: get_adj_facets_on_shell (PRIVATE)
1395 : : //Description: non recursive function that creates a list of all facets connected
1396 : : // the passed in face that are on the same surface
1397 : : //Author: sjowen
1398 : : //Date: 12/22/00
1399 : : //=============================================================================
1400 : 11 : CubitStatus FacetDataUtil::get_adj_facets_on_shell(
1401 : : CubitFacet *start_facet_ptr,
1402 : : DLIList<CubitFacet *> *shell_ptr,
1403 : : CubitBoolean &is_water_tight,
1404 : : int mydebug)
1405 : : {
1406 : 11 : int found = 0;
1407 : : int ii, jj;
1408 : 11 : CubitStatus stat = CUBIT_SUCCESS;
1409 [ + - ]: 11 : DLIList<CubitFacet*> temp_list;
1410 : 11 : CubitFacet *facet_ptr = NULL;
1411 : 11 : CubitFacet *adj_facet_ptr = NULL;
1412 [ + - ][ + - ]: 22 : DLIList<CubitFacetEdge *>edge_list;
1413 : 11 : CubitFacetEdge *edge_ptr = NULL;
1414 [ + - ][ + - ]: 22 : DLIList<CubitFacet *>adj_facet_list;
1415 : :
1416 : : // add this facet to the list
1417 : :
1418 [ + - ]: 11 : temp_list.append( start_facet_ptr );
1419 [ + - ]: 11 : start_facet_ptr->marked( 0 );
1420 : :
1421 [ + - ][ + + ]: 143 : while (temp_list.size())
1422 : : {
1423 [ + - ]: 132 : facet_ptr = temp_list.pop();
1424 [ + - ][ + - ]: 132 : if (facet_ptr->marked() == 0)
1425 : : {
1426 [ + - ]: 132 : shell_ptr->append( facet_ptr );
1427 [ - + ]: 132 : if (mydebug)
1428 : : {
1429 [ # # ]: 0 : GfxDebug::draw_facet(facet_ptr, CUBIT_RED_INDEX);
1430 [ # # ]: 0 : GfxDebug::flush();
1431 : : }
1432 [ + - ]: 132 : edge_list.clean_out();
1433 [ + - ]: 132 : facet_ptr->edges( edge_list );
1434 [ + - ][ + + ]: 528 : for (ii=0; ii<edge_list.size(); ii++)
1435 : : {
1436 [ + - ]: 396 : edge_ptr = edge_list.get_and_step();
1437 [ + - ]: 396 : adj_facet_list.clean_out();
1438 [ + - ]: 396 : edge_ptr->facets( adj_facet_list );
1439 : 396 : found = 0;
1440 : :
1441 [ + - ][ + + ]: 1155 : for (jj=0; jj<adj_facet_list.size() && !found; jj++)
[ + + ][ + + ]
1442 : : {
1443 [ + - ]: 759 : adj_facet_ptr = adj_facet_list.get_and_step();
1444 [ + + ]: 759 : if (adj_facet_ptr != facet_ptr)
1445 : : {
1446 : :
1447 : : // go to its neighbor if it is part of the surface
1448 : :
1449 [ + - ][ + + ]: 396 : if (adj_facet_ptr->marked() == 1)
1450 : : {
1451 [ + - ]: 121 : temp_list.append( adj_facet_ptr );
1452 [ + - ]: 121 : adj_facet_ptr->marked( 0 );
1453 : 121 : found = 1;
1454 : : }
1455 : : }
1456 : :
1457 : : }
1458 [ + - ][ + - ]: 396 : if (is_water_tight && adj_facet_list.size() == 1)
[ - + ][ - + ]
1459 : : {
1460 : 0 : is_water_tight = CUBIT_FALSE;
1461 : : }
1462 : : }
1463 : : }
1464 : : }
1465 : :
1466 [ + - ]: 11 : return stat;
1467 : : }
1468 : :
1469 : : //=============================================================================
1470 : : //Function: stitch_facets (PRIVATE)
1471 : : //Description: attempt to merge facets to form a watertight model
1472 : : //Author: sjowen
1473 : : //Date: 9/9/03
1474 : : //=============================================================================
1475 : 0 : CubitStatus FacetDataUtil::stitch_facets(
1476 : : DLIList<DLIList<CubitFacet *> *> &shell_list,
1477 : : double tol,
1478 : : CubitBoolean &is_water_tight,
1479 : : CubitBoolean write_result)
1480 : : {
1481 : :
1482 : 0 : CubitStatus rv = CUBIT_SUCCESS;
1483 : 0 : int npmerge = 0;
1484 : 0 : int nemerge = 0;
1485 [ # # ]: 0 : DLIList<CubitPoint*> unmerged_points;
1486 : 0 : is_water_tight = CUBIT_FALSE;
1487 : :
1488 [ # # ]: 0 : rv = merge_coincident_vertices( shell_list, tol, npmerge, nemerge, unmerged_points );
1489 [ # # ]: 0 : if (rv != CUBIT_SUCCESS)
1490 : 0 : return rv;
1491 : :
1492 [ # # ]: 0 : int nnomerge = unmerged_points.size();
1493 [ # # ]: 0 : if (nnomerge == 0)
1494 : : {
1495 : 0 : is_water_tight = CUBIT_TRUE;
1496 : : }
1497 : : else
1498 : : {
1499 : : rv = merge_colinear_vertices( shell_list, tol, unmerged_points,
1500 [ # # ]: 0 : npmerge, nemerge, nnomerge );
1501 [ # # ]: 0 : if (rv != CUBIT_SUCCESS)
1502 : 0 : return rv;
1503 : :
1504 [ # # ]: 0 : if (nnomerge == 0)
1505 : : {
1506 : 0 : is_water_tight = CUBIT_TRUE;
1507 : 0 : int mydebug = 0;
1508 [ # # ]: 0 : if (mydebug) // make sure its really water-tight
1509 : : {
1510 : : DLIList<CubitFacet *> *shell_ptr;
1511 [ # # ]: 0 : DLIList<CubitFacetEdge *>boundary_edges;
1512 [ # # ][ # # ]: 0 : for (int ii=0; ii<shell_list.size(); ii++)
1513 : : {
1514 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
1515 [ # # ]: 0 : FacetDataUtil::get_boundary_edges(*shell_ptr, boundary_edges);
1516 : : }
1517 [ # # ][ # # ]: 0 : if (boundary_edges.size() > 0)
1518 : : {
1519 [ # # ][ # # ]: 0 : PRINT_ERROR("Not Water-tight!\n");
[ # # ][ # # ]
1520 [ # # ]: 0 : }
1521 : : }
1522 : : }
1523 : : }
1524 : :
1525 [ # # ]: 0 : if (write_result)
1526 : : {
1527 [ # # ][ # # ]: 0 : if (npmerge > 0 || nemerge > 0)
1528 : : {
1529 [ # # ][ # # ]: 0 : PRINT_INFO("%d facet vertices and %d facet edges were successfully merged.\n",
[ # # ]
1530 [ # # ]: 0 : npmerge, nemerge);
1531 : : }
1532 [ # # ]: 0 : if (is_water_tight)
1533 : : {
1534 [ # # ][ # # ]: 0 : PRINT_INFO("Facets are water-tight.\n");
[ # # ][ # # ]
1535 : : }
1536 : : else
1537 : : {
1538 [ # # ][ # # ]: 0 : PRINT_INFO("%d facet vertices on model boundary detected that could not be merged.\n", nnomerge);
[ # # ][ # # ]
1539 : : }
1540 : : }
1541 : :
1542 [ # # ]: 0 : return rv;
1543 : : }
1544 : :
1545 : : //=============================================================================
1546 : : //Function: merge_colinear_vertices (PRIVATE)
1547 : : //Description: check if any vertices fall on a free facet edge. If so - split
1548 : : // the adjacent facet to merge
1549 : : //Author: sjowen
1550 : : //Date: 1/19/2004
1551 : : //=============================================================================
1552 : 0 : CubitStatus FacetDataUtil::merge_colinear_vertices(
1553 : : DLIList<DLIList<CubitFacet *> *> &shell_list,
1554 : : double tol,
1555 : : DLIList<CubitPoint *> &merge_point_list, // points to attempt to merge
1556 : : int &npmerge, // return number of vertices merged
1557 : : int &nemerge, // return number of edges merged
1558 : : int &nnomerge) // return number of vertices in list NOT merged
1559 : : {
1560 : 0 : nnomerge = 0;
1561 : 0 : int mydebug = 0;
1562 : :
1563 : : int ii, jj, kk, ll, mm, pt_shell_id, edge_shell_id;
1564 [ # # ]: 0 : RTree <CubitFacetEdge*> r_tree(tol);
1565 [ # # ]: 0 : shell_list.reset();
1566 : 0 : DLIList<CubitFacet *> *shell_ptr = NULL;
1567 : : CubitPoint *pt;
1568 : : CubitFacetEdge *edge;
1569 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *>boundary_edges;
1570 [ # # ][ # # ]: 0 : for (ii=0; ii<shell_list.size(); ii++)
1571 : : {
1572 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
1573 : :
1574 : : // get the boundary edges from the shell - these are candidates for merging.
1575 : : // note that this won't merge internal facets
1576 : :
1577 [ # # ]: 0 : DLIList<CubitFacetEdge *>shell_boundary_edges;
1578 [ # # ]: 0 : FacetDataUtil::get_boundary_edges(*shell_ptr, shell_boundary_edges);
1579 : :
1580 : : // mark each of the edges with a shell id so we know when to merge shells
1581 : : // and add them to the kdtree for fast spatial searching
1582 : 0 : edge_shell_id = ii+1;
1583 [ # # ][ # # ]: 0 : for(jj=0; jj<shell_boundary_edges.size(); jj++)
1584 : : {
1585 [ # # ]: 0 : edge = shell_boundary_edges.get_and_step();
1586 [ # # ][ # # ]: 0 : edge->point(0)->marked(edge_shell_id);
1587 [ # # ][ # # ]: 0 : edge->point(1)->marked(edge_shell_id);
1588 [ # # ]: 0 : edge->marked(edge_shell_id);
1589 [ # # ]: 0 : r_tree.add(edge);
1590 [ # # ]: 0 : boundary_edges.append(edge);
1591 : : }
1592 [ # # ][ # # ]: 0 : for(jj=0; jj<shell_ptr->size(); jj++)
1593 [ # # ][ # # ]: 0 : shell_ptr->get_and_step()->marked( edge_shell_id );
1594 [ # # ]: 0 : }
1595 : :
1596 : : // find points in merge_point_list that are colinear with edges in boundary_edges
1597 : :
1598 [ # # ][ # # ]: 0 : CubitBox ptbox;
1599 [ # # ][ # # ]: 0 : CubitVector ptmin, ptmax, coord;
[ # # ]
1600 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *>close_edges;
1601 : : CubitFacetEdge *close_edge;
1602 : : CubitFacet *facet;
1603 : : CubitPoint *p0, *p1;
1604 [ # # ][ # # ]: 0 : DLIList<CubitPoint*>adj_pt_list;
1605 [ # # ][ # # ]: 0 : DLIList<CubitPoint*>del_points;
1606 [ # # ][ # # ]: 0 : for(ii=0; ii<merge_point_list.size(); ii++)
1607 : : {
1608 [ # # ]: 0 : pt = merge_point_list.get_and_step();
1609 [ # # ][ # # ]: 0 : if (pt->marked() < 0) // has already been merged
1610 : 0 : continue;
1611 : :
1612 : : // find the closest edges in the rtree
1613 : :
1614 [ # # ][ # # ]: 0 : coord = pt->coordinates();
1615 [ # # ][ # # ]: 0 : ptmin.set( coord.x() - tol, coord.y() - tol, coord.z() - tol );
[ # # ][ # # ]
1616 [ # # ][ # # ]: 0 : ptmax.set( coord.x() + tol, coord.y() + tol, coord.z() + tol );
[ # # ][ # # ]
1617 [ # # ]: 0 : ptbox.reset( ptmin, ptmax );
1618 [ # # ]: 0 : close_edges.clean_out();
1619 [ # # ]: 0 : r_tree.find(ptbox, close_edges);
1620 : :
1621 : :
1622 : : //remove any close edges that already share the merge point
1623 [ # # ]: 0 : DLIList<CubitFacetEdge*> adj_edges;
1624 [ # # ][ # # ]: 0 : for (jj=close_edges.size(); jj--; )
1625 : : {
1626 [ # # ]: 0 : close_edge = close_edges.get();
1627 : :
1628 [ # # ][ # # ]: 0 : if (close_edge->point(0) == pt || close_edge->point(1) == pt)
[ # # ][ # # ]
[ # # ]
1629 : : {
1630 [ # # ]: 0 : close_edges.change_to( NULL );
1631 [ # # ]: 0 : adj_edges.append( close_edge );
1632 : : }
1633 : :
1634 [ # # ]: 0 : close_edges.step();
1635 : : }
1636 [ # # ]: 0 : close_edges.remove_all_with_value( NULL );
1637 : :
1638 [ # # ][ # # ]: 0 : if( close_edges.size() == 0 && adj_edges.size() > 0 )
[ # # ][ # # ]
[ # # ]
1639 : : {
1640 : : //use a coareser tolerance to get some edges
1641 : : //one tenth the average length of attached edges
1642 : 0 : double temp_tol = 0;
1643 [ # # ][ # # ]: 0 : for( int i=adj_edges.size(); i--; )
1644 [ # # ][ # # ]: 0 : temp_tol += adj_edges.get_and_step()->length();
1645 : :
1646 [ # # ]: 0 : temp_tol /= adj_edges.size();
1647 : 0 : temp_tol *= 0.1;
1648 : :
1649 : : //bump up tolerance to get more edges
1650 [ # # ][ # # ]: 0 : ptmin.set( coord.x() - temp_tol, coord.y() - temp_tol, coord.z() - temp_tol );
[ # # ][ # # ]
1651 [ # # ][ # # ]: 0 : ptmax.set( coord.x() + temp_tol, coord.y() + temp_tol, coord.z() + temp_tol );
[ # # ][ # # ]
1652 [ # # ]: 0 : ptbox.reset( ptmin, ptmax );
1653 [ # # ]: 0 : close_edges.clean_out();
1654 [ # # ]: 0 : r_tree.find(ptbox, close_edges);
1655 : : }
1656 : :
1657 : :
1658 : : // We did find something - go try to merge
1659 : :
1660 : 0 : CubitBoolean was_merged = CUBIT_FALSE;
1661 [ # # ][ # # ]: 0 : for (jj=0; jj<close_edges.size() && !was_merged; jj++)
[ # # ][ # # ]
1662 : : {
1663 : : // test the next edge on the list
1664 : :
1665 [ # # ]: 0 : close_edge = close_edges.get_and_step();
1666 : :
1667 : : // make sure the edge does not contain the merge point
1668 : :
1669 [ # # ][ # # ]: 0 : if (close_edge->point(0) == pt || close_edge->point(1) == pt)
[ # # ][ # # ]
[ # # ]
1670 : 0 : continue;
1671 : :
1672 : : // check to see if the edge is within tolerance of the merge point
1673 : :
1674 [ # # ]: 0 : CubitVector loc_on_edge;
1675 : : CubitBoolean is_outside_edge;
1676 : 0 : double dist = close_edge->dist_to_edge(pt->coordinates(),
1677 [ # # ][ # # ]: 0 : loc_on_edge, is_outside_edge);
1678 : :
1679 : : // allow for some surface curvature. permit the point to be close to
1680 : : // the edge but not on. May want to modify this factor if it not closing
1681 : : // the facets correctly. If it's too big, it may get edges that aren't
1682 : : // really adjacent.
1683 [ # # ]: 0 : double edge_tol = 0.1 * close_edge->length();
1684 [ # # ][ # # ]: 0 : if (is_outside_edge || dist > edge_tol)
1685 : 0 : continue;
1686 : :
1687 : : // go merge the point with the edge
1688 : 0 : was_merged = CUBIT_TRUE;
1689 : :
1690 [ # # ]: 0 : if (mydebug)
1691 : : {
1692 [ # # ]: 0 : DLIList<CubitFacet *> fl;
1693 [ # # ]: 0 : pt->facets( fl );
1694 [ # # ]: 0 : dcolor(CUBIT_YELLOW_INDEX);
1695 [ # # ]: 0 : dfldraw(fl);
1696 [ # # ]: 0 : dcolor(CUBIT_BLUE_INDEX);
1697 [ # # ][ # # ]: 0 : dpoint(pt->coordinates());
1698 [ # # ]: 0 : CubitFacet *f = close_edge->adj_facet(0);
1699 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
1700 [ # # ]: 0 : dfdraw(f);
1701 [ # # ][ # # ]: 0 : dview();
1702 : : }
1703 : :
1704 : : // remove close_edge from the rtree
1705 : :
1706 [ # # ]: 0 : r_tree.remove( close_edge );
1707 [ # # ]: 0 : edge_shell_id = close_edge->marked();
1708 [ # # ]: 0 : close_edge->marked(0);
1709 : :
1710 : : // split the edge
1711 : :
1712 [ # # ][ # # ]: 0 : assert(1 == close_edge->num_adj_facets()); // assumes we are splitting a boundary facet
1713 [ # # ]: 0 : facet = close_edge->adj_facet(0);
1714 [ # # ]: 0 : CubitFacetData *dfacet = dynamic_cast<CubitFacetData *>(facet);
1715 : 0 : CubitPoint *new_pt = dfacet->split_edge(close_edge->point(0),
1716 : 0 : close_edge->point(1),
1717 [ # # ][ # # ]: 0 : pt->coordinates());
[ # # ][ # # ]
1718 [ # # ]: 0 : int facet_tool_id = facet->tool_id();
1719 [ # # ]: 0 : new_pt->marked(edge_shell_id);
1720 : :
1721 : : // add any new facets to the shell and create missing edges
1722 : :
1723 [ # # ]: 0 : for (ll=0; ll<edge_shell_id; ll++)
1724 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
1725 [ # # ]: 0 : DLIList<CubitFacet *>adj_facets;
1726 [ # # ]: 0 : new_pt->facets( adj_facets );
1727 [ # # ][ # # ]: 0 : for (ll=0; ll<adj_facets.size(); ll++)
1728 : : {
1729 [ # # ]: 0 : facet = adj_facets.get_and_step();
1730 [ # # ][ # # ]: 0 : if (!facet->marked())
1731 : : {
1732 [ # # ]: 0 : facet->marked(edge_shell_id);
1733 [ # # ]: 0 : shell_ptr->append(facet);
1734 [ # # ]: 0 : for (mm=0; mm<3; mm++) {
1735 [ # # ][ # # ]: 0 : if (!(edge = facet->edge(mm)))
1736 : : {
1737 [ # # ]: 0 : facet->get_edge_pts(mm, p0, p1);
1738 [ # # ][ # # ]: 0 : edge = (CubitFacetEdge *) new CubitFacetEdgeData( p0, p1 );
1739 [ # # ]: 0 : edge->marked( 0 );
1740 : : }
1741 : : }
1742 [ # # ]: 0 : facet->set_tool_id(facet_tool_id);
1743 : : }
1744 : : }
1745 : :
1746 : : // merge the points,
1747 : :
1748 [ # # ]: 0 : merge_points( pt, new_pt, nemerge, &r_tree );
1749 : 0 : npmerge++;
1750 : :
1751 : : // add any new edges to the rtree
1752 : :
1753 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *> adj_edges;
[ # # ]
1754 [ # # ]: 0 : pt->edges( adj_edges );
1755 [ # # ][ # # ]: 0 : for (kk=0; kk<adj_edges.size(); kk++)
1756 : : {
1757 [ # # ]: 0 : CubitFacetEdge *adj_edge = adj_edges.get_and_step();
1758 [ # # ][ # # ]: 0 : if (!adj_edge->marked() && adj_edge->num_adj_facets() == 1)
[ # # ][ # # ]
[ # # ]
1759 : : {
1760 [ # # ][ # # ]: 0 : adj_edge->marked(pt->marked());
1761 [ # # ]: 0 : r_tree.add(adj_edge);
1762 : : }
1763 : : }
1764 : :
1765 : : // see if shells need to merge and then merge
1766 : :
1767 [ # # ]: 0 : pt_shell_id = pt->marked();
1768 [ # # ]: 0 : if (pt_shell_id != edge_shell_id)
1769 : : {
1770 : :
1771 : : // get the shell containing the close point. Nullify the
1772 : : // pointer in the shell list and move all of its facets
1773 : : // to the other shell.
1774 : :
1775 : 0 : int delete_shell = edge_shell_id;
1776 : 0 : DLIList<CubitFacet *> *delete_shell_ptr = NULL;
1777 [ # # ]: 0 : shell_list.reset();
1778 [ # # ]: 0 : for (ll=0; ll<delete_shell; ll++)
1779 [ # # ]: 0 : delete_shell_ptr = shell_list.get_and_step();
1780 [ # # ]: 0 : shell_list.back();
1781 [ # # ]: 0 : shell_list.change_to(NULL);
1782 : :
1783 [ # # ]: 0 : if(!delete_shell_ptr)
1784 : 0 : return CUBIT_FAILURE;
1785 : :
1786 : : // mark all the points on the delete_shell as now being part of
1787 : : // the other shell.
1788 : :
1789 [ # # ][ # # ]: 0 : for(ll=0; ll<delete_shell_ptr->size(); ll++)
1790 : : {
1791 [ # # ]: 0 : facet = delete_shell_ptr->get_and_step();
1792 [ # # ]: 0 : facet->marked( pt_shell_id );
1793 [ # # ][ # # ]: 0 : facet->point( 0 )->marked( pt_shell_id );
1794 [ # # ][ # # ]: 0 : facet->point( 1 )->marked( pt_shell_id );
1795 [ # # ][ # # ]: 0 : facet->point( 2 )->marked( pt_shell_id );
1796 [ # # ][ # # ]: 0 : facet->edge( 0 )->marked( pt_shell_id );
1797 [ # # ][ # # ]: 0 : facet->edge( 1 )->marked( pt_shell_id );
1798 [ # # ][ # # ]: 0 : facet->edge( 2 )->marked( pt_shell_id );
1799 : : }
1800 : :
1801 : : // append the two lists together
1802 : :
1803 [ # # ]: 0 : shell_list.reset();
1804 [ # # ]: 0 : for (ll=0; ll<pt_shell_id; ll++)
1805 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
1806 [ # # ]: 0 : *shell_ptr += (*delete_shell_ptr);
1807 [ # # ][ # # ]: 0 : delete delete_shell_ptr;
1808 : : }
1809 : :
1810 : : // set the marked flag to negative to indicate that it has been
1811 : : // merged and it needs to be deleted.
1812 : :
1813 [ # # ]: 0 : del_points.append(new_pt);
1814 [ # # ][ # # ]: 0 : new_pt->marked( -new_pt->marked() );
1815 [ # # ][ # # ]: 0 : pt->marked( -pt->marked() );
[ # # ][ # # ]
1816 : 0 : } // close_edges
1817 : :
1818 [ # # ]: 0 : if (!was_merged)
1819 : : {
1820 : 0 : nnomerge++;
1821 [ # # ]: 0 : if (mydebug)
1822 : : {
1823 [ # # ]: 0 : dcolor(CUBIT_BLUE_INDEX);
1824 [ # # ]: 0 : dpdraw(pt);
1825 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
1826 [ # # ]: 0 : CubitVector mmin(-1e10, -1e10, -1e10);
1827 [ # # ]: 0 : CubitVector mmax(1e10, 1e10, 1e10);
1828 [ # # ]: 0 : ptbox.reset(mmin, mmax);
1829 [ # # ]: 0 : r_tree.find(ptbox, close_edges);
1830 [ # # ][ # # ]: 0 : for(int zz=0; zz<close_edges.size(); zz++)
1831 : : {
1832 [ # # ]: 0 : edge = close_edges.get_and_step();
1833 [ # # ]: 0 : dedraw(edge);
1834 [ # # ]: 0 : CubitVector loc_on_edge;
1835 : : CubitBoolean is_on_edge;
1836 [ # # ][ # # ]: 0 : double dist = edge->dist_to_edge(pt->coordinates(), loc_on_edge, is_on_edge);
1837 [ # # ][ # # ]: 0 : PRINT_INFO("edge %d dist = %f is_on_edge = %d\n", edge->id(), dist, is_on_edge);
[ # # ][ # # ]
[ # # ]
1838 : : }
1839 [ # # ][ # # ]: 0 : dview();
[ # # ]
1840 : : }
1841 : : }
1842 : 0 : } // merge_point_list
1843 : :
1844 : : // compress the shell list
1845 : :
1846 [ # # ]: 0 : shell_list.remove_all_with_value(NULL);
1847 : :
1848 : : // delete points that were merged
1849 : :
1850 [ # # ][ # # ]: 0 : for (ii=0; ii<del_points.size(); ii++)
1851 : : {
1852 [ # # ]: 0 : pt = del_points.get_and_step();
1853 [ # # ][ # # ]: 0 : if (pt->marked() < 0)
1854 : : {
1855 [ # # ][ # # ]: 0 : assert( pt->num_adj_facets() == 0 );
1856 [ # # ][ # # ]: 0 : delete pt;
1857 : : }
1858 : : }
1859 : :
1860 [ # # ]: 0 : return CUBIT_SUCCESS;
1861 : : }
1862 : :
1863 : : //=============================================================================
1864 : : //Function: merge_points (PRIVATE)
1865 : : //Description: merge two cubit points where it has been determined that they
1866 : : // are coincident. Also handle merging their adjacent edges where
1867 : : // appropriate
1868 : : //Author: sjowen
1869 : : //Date: 1/19/2004
1870 : : //=============================================================================
1871 : 0 : CubitStatus FacetDataUtil::merge_points(
1872 : : CubitPoint *pt0, CubitPoint *pt1,
1873 : : int &nemerge, RTree <CubitFacetEdge*> *r_tree) // number of edges we had to merge top do this
1874 : : {
1875 : 0 : int mydebug = 0;
1876 : :
1877 : : // merge the points
1878 : 0 : pt0->merge_points( pt1, CUBIT_TRUE );
1879 : : //pt0->set_as_feature();
1880 : :
1881 [ # # ]: 0 : if (mydebug)
1882 : : {
1883 [ # # ]: 0 : DLIList<CubitFacet *>adj_facets;
1884 [ # # ]: 0 : pt0->facets( adj_facets );
1885 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
1886 [ # # ][ # # ]: 0 : dfldraw(adj_facets);
1887 : : }
1888 : :
1889 : : // merge edges
1890 : :
1891 : : int ll, mm;
1892 : : bool edge_merged;
1893 [ # # ]: 0 : do {
1894 : 0 : edge_merged = false;
1895 : : CubitFacetEdge *edge, *other_edge;
1896 : :
1897 [ # # ]: 0 : DLIList<CubitFacetEdge *> adj_edges;
1898 [ # # ]: 0 : pt0->edges( adj_edges );
1899 [ # # ][ # # ]: 0 : for(ll=0; ll<adj_edges.size() && !edge_merged; ll++)
[ # # ][ # # ]
1900 : : {
1901 [ # # ]: 0 : edge = adj_edges.get_and_step();
1902 [ # # ][ # # ]: 0 : for(mm=0; mm<adj_edges.size() && !edge_merged; mm++)
[ # # ][ # # ]
1903 : : {
1904 [ # # ]: 0 : other_edge = adj_edges.get_and_step();
1905 [ # # ][ # # ]: 0 : if (other_edge != edge &&
[ # # ]
1906 [ # # ][ # # ]: 0 : edge->other_point(pt0) == other_edge->other_point(pt0))
1907 : : {
1908 [ # # ]: 0 : CubitFacetEdgeData *dedge = dynamic_cast<CubitFacetEdgeData*>(edge);
1909 [ # # ]: 0 : CubitFacetEdgeData *dother_edge = dynamic_cast<CubitFacetEdgeData*>(other_edge);
1910 : : // mbrewer (11/16/2005 for Bug 5049)
1911 [ # # ]: 0 : if(r_tree){
1912 : : // r_tree->remove(dother_edge);
1913 [ # # ]: 0 : CubitBoolean temp_bool = r_tree->remove(dother_edge);
1914 [ # # ]: 0 : if(!temp_bool){
1915 [ # # ][ # # ]: 0 : PRINT_DEBUG_139("FacetDataUtil: D_OTHER_EDGE did not exist in RTREE.\n");
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1916 : : }
1917 : : }
1918 [ # # ][ # # ]: 0 : if (dedge->merge_edges( dother_edge ) == CUBIT_SUCCESS)
1919 : : {
1920 : 0 : nemerge++;
1921 : 0 : edge_merged = true;
1922 [ # # ]: 0 : dedge->set_as_feature();
1923 : : }
1924 : : }
1925 : : }
1926 [ # # ]: 0 : }
1927 : : } while (edge_merged);
1928 : 0 : return CUBIT_SUCCESS;
1929 : : }
1930 : :
1931 : :
1932 : 0 : CubitStatus FacetDataUtil::merge_coincident_vertices( DLIList<CubitPoint*> &points )
1933 : : {
1934 : 0 : points.reset();
1935 : 0 : CubitPoint *surviving_point = points.pop();
1936 [ # # ]: 0 : while( points.size() )
1937 : : {
1938 : 0 : int nemerge = 0;
1939 [ # # ][ # # ]: 0 : merge_points( surviving_point, points.pop(), nemerge );
1940 : : }
1941 : :
1942 : 0 : return CUBIT_SUCCESS;
1943 : : }
1944 : :
1945 : :
1946 : :
1947 : : //=============================================================================
1948 : : //Function: merge_coincident_vertices (PRIVATE)
1949 : : //Description: merge vertices (and connected facets and edges) if they are
1950 : : // within tolerance.
1951 : : //Author: sjowen
1952 : : //Date: 9/9/03
1953 : : //=============================================================================
1954 : 0 : CubitStatus FacetDataUtil::merge_coincident_vertices(
1955 : : DLIList<DLIList<CubitFacet *> *> &shell_list,
1956 : : double tol,
1957 : : int &npmerge, // return number of vertices merged
1958 : : int &nemerge, // return number of edges merged
1959 : : DLIList<CubitPoint *> &unmerged_points) // return the vertices on boundary NOT merged
1960 : : {
1961 : 0 : int mydebug = 0;
1962 : 0 : npmerge = 0;
1963 : 0 : nemerge = 0;
1964 : : int ii, jj, kk, ll, shell_id;
1965 [ # # ]: 0 : KDDTree <CubitPoint*> kd_tree(tol, false);
1966 [ # # ]: 0 : CubitPoint::set_box_tol( tol );
1967 [ # # ]: 0 : shell_list.reset();
1968 : 0 : DLIList<CubitFacet *> *shell_ptr = NULL;
1969 : : CubitPoint *pt;
1970 [ # # ][ # # ]: 0 : DLIList<CubitPoint *>boundary_points;
1971 [ # # ][ # # ]: 0 : DLIList<CubitPoint *> del_points;
1972 [ # # ][ # # ]: 0 : for (ii=0; ii<shell_list.size(); ii++)
1973 : : {
1974 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
1975 : :
1976 : : // get the boundary points from the shell - these are candidates for merging.
1977 : : // note that this won't merge internal facets
1978 : :
1979 [ # # ]: 0 : DLIList<CubitPoint *>shell_boundary_points;
1980 [ # # ]: 0 : FacetDataUtil::get_boundary_points(*shell_ptr, shell_boundary_points);
1981 : :
1982 : : // mark each of the points with a shell id so we know when to merge shells
1983 : : // and add them to the kdtree for fast spatial searching
1984 : :
1985 : 0 : shell_id = ii+1;
1986 [ # # ][ # # ]: 0 : for(jj=0; jj<shell_boundary_points.size(); jj++)
1987 : : {
1988 [ # # ]: 0 : pt = shell_boundary_points.get_and_step();
1989 [ # # ]: 0 : pt->marked(shell_id);
1990 [ # # ]: 0 : kd_tree.add(pt);
1991 [ # # ]: 0 : boundary_points.append(pt);
1992 : : }
1993 [ # # ]: 0 : }
1994 [ # # ]: 0 : kd_tree.balance();
1995 : :
1996 : : // find points that are coincident
1997 : :
1998 [ # # ][ # # ]: 0 : CubitBox ptbox;
1999 [ # # ][ # # ]: 0 : CubitVector ptmin, ptmax, coord;
[ # # ]
2000 [ # # ][ # # ]: 0 : DLIList<CubitPoint *>close_points;
2001 : 0 : CubitPoint *close_pt = NULL;
2002 : : CubitFacet *facet;
2003 [ # # ][ # # ]: 0 : DLIList<CubitPoint*>adj_pt_list;
2004 : :
2005 [ # # ][ # # ]: 0 : for(ii=0; ii<boundary_points.size(); ii++)
2006 : : {
2007 [ # # ]: 0 : pt = boundary_points.get_and_step();
2008 [ # # ][ # # ]: 0 : if (pt->marked() < 0) // has already been merged
2009 : 0 : continue;
2010 : :
2011 : : // find the closest points in the kdtree
2012 : :
2013 [ # # ][ # # ]: 0 : coord = pt->coordinates();
2014 [ # # ][ # # ]: 0 : ptmin.set( coord.x() - tol, coord.y() - tol, coord.z() - tol );
[ # # ][ # # ]
2015 [ # # ][ # # ]: 0 : ptmax.set( coord.x() + tol, coord.y() + tol, coord.z() + tol );
[ # # ][ # # ]
2016 [ # # ]: 0 : ptbox.reset( ptmin, ptmax );
2017 [ # # ]: 0 : close_points.clean_out();
2018 [ # # ]: 0 : kd_tree.find(ptbox, close_points);
2019 : :
2020 : : // if it didn't find anything to merge with, then we aren't water-tight
2021 : :
2022 : 0 : CubitBoolean was_merged = CUBIT_FALSE;
2023 : :
2024 : : // We did find something - go try to merge
2025 : :
2026 [ # # ][ # # ]: 0 : for (jj=0; jj<close_points.size(); jj++)
2027 : : {
2028 [ # # ]: 0 : close_pt = close_points.get_and_step();
2029 [ # # ]: 0 : if (close_pt == pt)
2030 : 0 : continue;
2031 [ # # ][ # # ]: 0 : if (close_pt->marked() < 0) // has already been merged
2032 : 0 : continue;
2033 [ # # ][ # # ]: 0 : assert(close_points.size() >= 1);
2034 : :
2035 : : // make sure this point is not already one of its neighbors
2036 : : // so we don't collapse a triangle
2037 : :
2038 : 0 : CubitBoolean is_adjacent = CUBIT_FALSE;
2039 [ # # ]: 0 : adj_pt_list.clean_out();
2040 [ # # ]: 0 : close_pt->adjacent_points( adj_pt_list );
2041 [ # # ][ # # ]: 0 : for (kk=0; kk<adj_pt_list.size() && !is_adjacent; kk++)
[ # # ][ # # ]
2042 : : {
2043 [ # # ][ # # ]: 0 : if (adj_pt_list.get_and_step() == pt)
2044 : 0 : is_adjacent = CUBIT_TRUE;
2045 : : }
2046 [ # # ]: 0 : if (!is_adjacent)
2047 : : {
2048 : : // merge the points
2049 : :
2050 [ # # ]: 0 : merge_points( pt, close_pt, nemerge );
2051 : 0 : npmerge++;
2052 : 0 : was_merged = CUBIT_TRUE;
2053 : :
2054 : : // see if shells need to merge and then merge
2055 : :
2056 [ # # ]: 0 : shell_id = pt->marked();
2057 [ # # ][ # # ]: 0 : if (shell_id != close_pt->marked())
2058 : : {
2059 : :
2060 : : // get the shell containing the close point. Nullify the
2061 : : // pointer in the shell list and move all of its facets
2062 : : // to the other shell.
2063 : :
2064 [ # # ]: 0 : int delete_shell = close_pt->marked();
2065 : 0 : DLIList<CubitFacet *> *delete_shell_ptr = NULL;
2066 [ # # ]: 0 : shell_list.reset();
2067 [ # # ]: 0 : for (ll=0; ll<delete_shell; ll++)
2068 [ # # ]: 0 : delete_shell_ptr = shell_list.get_and_step();
2069 [ # # ]: 0 : shell_list.back();
2070 [ # # ]: 0 : shell_list.change_to(NULL);
2071 : :
2072 : : // mark all the points on the delete_shell as now being part of
2073 : : // the other shell.
2074 : :
2075 [ # # ][ # # ]: 0 : for(ll=0; ll<delete_shell_ptr->size(); ll++)
2076 : : {
2077 [ # # ]: 0 : facet = delete_shell_ptr->get_and_step();
2078 [ # # ][ # # ]: 0 : facet->point( 0 )->marked( shell_id );
2079 [ # # ][ # # ]: 0 : facet->point( 1 )->marked( shell_id );
2080 [ # # ][ # # ]: 0 : facet->point( 2 )->marked( shell_id );
2081 : : }
2082 : :
2083 : : // append the two lists together
2084 : :
2085 [ # # ]: 0 : shell_list.reset();
2086 [ # # ]: 0 : for (ll=0; ll<shell_id; ll++)
2087 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
2088 [ # # ]: 0 : *shell_ptr += (*delete_shell_ptr);
2089 [ # # ][ # # ]: 0 : delete delete_shell_ptr;
2090 : : }
2091 : :
2092 : : // set the marked flag to negative to indicate that it has been
2093 : : // merged and it need to be deleted.
2094 : :
2095 [ # # ][ # # ]: 0 : if (close_pt->marked() > 0)
2096 [ # # ][ # # ]: 0 : close_pt->marked( -close_pt->marked() );
2097 [ # # ]: 0 : del_points.append( close_pt );
2098 : : }
2099 : : }
2100 [ # # ]: 0 : if (was_merged)
2101 : : {
2102 [ # # ][ # # ]: 0 : if (pt->marked() > 0)
2103 [ # # ][ # # ]: 0 : pt->marked( -pt->marked() );
2104 : : }
2105 : : else
2106 : : {
2107 : : // check to see if it was already merged
2108 [ # # ][ # # ]: 0 : if (close_points.size() == 1 && close_pt == pt)
[ # # ][ # # ]
2109 : : {
2110 : : CubitFacetEdge *edge_ptr;
2111 [ # # ]: 0 : DLIList<CubitFacetEdge *>adj_edges;
2112 [ # # ]: 0 : pt->edges(adj_edges);
2113 : 0 : CubitBoolean on_boundary = CUBIT_FALSE;
2114 [ # # ][ # # ]: 0 : for(kk=0; kk<adj_edges.size() && !on_boundary; kk++)
[ # # ][ # # ]
2115 : : {
2116 [ # # ]: 0 : edge_ptr = adj_edges.get_and_step();
2117 [ # # ][ # # ]: 0 : if (edge_ptr->num_adj_facets() == 1)
2118 : 0 : on_boundary = CUBIT_TRUE;
2119 : : }
2120 [ # # ]: 0 : if (on_boundary)
2121 : : {
2122 : : //PRINT_INFO("Merging 'boundary' points.\n");
2123 : : //if (pt->marked() > 0) pt->marked( -pt->marked() );
2124 [ # # ]: 0 : unmerged_points.append(pt);
2125 : : }
2126 : : else
2127 : : {
2128 [ # # ][ # # ]: 0 : if (pt->marked() > 0) pt->marked( -pt->marked() );
[ # # ][ # # ]
2129 [ # # ]: 0 : }
2130 : : }
2131 : : else
2132 : : {
2133 : : // otherwise save it as an unmerged point to be handled later
2134 [ # # ]: 0 : unmerged_points.append(pt);
2135 : : }
2136 : : }
2137 : : }
2138 : :
2139 : : // compress the shell list
2140 : :
2141 [ # # ]: 0 : shell_list.remove_all_with_value(NULL);
2142 : :
2143 : : // delete points that were merged
2144 : :
2145 [ # # ][ # # ]: 0 : for (ii=0; ii<del_points.size(); ii++)
2146 : : {
2147 [ # # ]: 0 : pt = del_points.get_and_step();
2148 [ # # ][ # # ]: 0 : if (pt->marked() < 0)
2149 : : {
2150 [ # # ][ # # ]: 0 : assert( pt->num_adj_facets() == 0 );
2151 [ # # ][ # # ]: 0 : delete pt;
2152 : : }
2153 : : }
2154 : :
2155 : : //double-check that these are still boundary points.
2156 : : //found case where two shells were tied together by a single facet point.
2157 : : //this point was initially deemed an unmerged boundary point but later
2158 : : //adjacent boundary edges got merged so it should not be a boundary
2159 : : //point anymore.
2160 [ # # ][ # # ]: 0 : for( int k=0; k<unmerged_points.size(); k++ )
2161 : : {
2162 [ # # ]: 0 : DLIList<CubitFacetEdge *>adj_edges;
2163 [ # # ][ # # ]: 0 : unmerged_points[k]->edges(adj_edges);
2164 : 0 : CubitBoolean on_boundary = CUBIT_FALSE;
2165 [ # # ][ # # ]: 0 : for(kk=0; kk<adj_edges.size(); kk++)
2166 : : {
2167 [ # # ]: 0 : CubitFacetEdge *edge_ptr = adj_edges.get_and_step();
2168 [ # # ][ # # ]: 0 : if (edge_ptr->num_adj_facets() == 1)
2169 : : {
2170 : 0 : on_boundary = CUBIT_TRUE;
2171 : 0 : break;
2172 : : }
2173 : : }
2174 [ # # ]: 0 : if( CUBIT_FALSE == on_boundary )
2175 [ # # ]: 0 : unmerged_points[k] = NULL;
2176 [ # # ]: 0 : }
2177 [ # # ]: 0 : unmerged_points.remove_all_with_value(NULL);
2178 : :
2179 [ # # ]: 0 : if (mydebug)
2180 : : {
2181 : : CubitFacetEdge *edge;
2182 [ # # ][ # # ]: 0 : for (ii=0; ii<shell_list.size(); ii++)
2183 : : {
2184 [ # # ]: 0 : shell_ptr = shell_list.get_and_step();
2185 [ # # ]: 0 : dcolor(CUBIT_GREEN_INDEX);
2186 [ # # ]: 0 : dfldraw(*shell_ptr);
2187 [ # # ]: 0 : dcolor(CUBIT_RED_INDEX);
2188 [ # # ][ # # ]: 0 : for(jj=0; jj<shell_ptr->size(); jj++)
2189 : : {
2190 [ # # ]: 0 : facet = shell_ptr->get_and_step();
2191 [ # # ]: 0 : for(kk=0; kk<3; kk++)
2192 : : {
2193 [ # # ]: 0 : edge = facet->edge(kk);
2194 [ # # ]: 0 : DLIList<FacetEntity *> myfacet_list;
2195 [ # # ]: 0 : edge->get_parents( myfacet_list );
2196 [ # # ][ # # ]: 0 : assert(myfacet_list.size() == edge->num_adj_facets());
[ # # ]
2197 [ # # ][ # # ]: 0 : if (myfacet_list.size() != 2)
2198 : : {
2199 [ # # ]: 0 : dedraw( edge );
2200 : : }
2201 [ # # ]: 0 : }
2202 : : }
2203 : : }
2204 [ # # ]: 0 : dcolor(CUBIT_BLUE_INDEX);
2205 [ # # ]: 0 : dpldraw(unmerged_points);
2206 [ # # ]: 0 : dview();
2207 : : }
2208 : :
2209 [ # # ]: 0 : return CUBIT_SUCCESS;
2210 : :
2211 : : }
2212 : :
2213 : : //=============================================================================
2214 : : //Function: delete_facets (PUBLIC)
2215 : : //Description: delete the facets and all associated edges and points from
2216 : : // a list of lists of facets (shell_list)
2217 : : //Author: sjowen
2218 : : //Date: 1/21/2004
2219 : : //=============================================================================
2220 : 0 : void FacetDataUtil::delete_facets(DLIList<DLIList<CubitFacet*>*> &shell_list)
2221 : : {
2222 : : int ii;
2223 [ # # ]: 0 : for (ii=0; ii<shell_list.size(); ii++)
2224 : : {
2225 : 0 : DLIList<CubitFacet*> *facet_list_ptr = shell_list.get_and_step();
2226 : 0 : delete_facets( *facet_list_ptr );
2227 : : }
2228 : 0 : }
2229 : :
2230 : : //=============================================================================
2231 : : //Function: delete_facets (PUBLIC)
2232 : : //Description: delete the facets and all associated edges and points from
2233 : : // a list of facets
2234 : : //Author: sjowen
2235 : : //Date: 1/21/2004
2236 : : //=============================================================================
2237 : 0 : void FacetDataUtil::delete_facets(DLIList<CubitFacet*> &facet_list)
2238 : : {
2239 : : int ii;
2240 : : CubitFacet *facet_ptr;
2241 [ # # ]: 0 : for (ii=0; ii<facet_list.size(); ii++)
2242 : : {
2243 : 0 : facet_ptr = facet_list.get_and_step();
2244 : 0 : delete_facet( facet_ptr );
2245 : : }
2246 : 0 : }
2247 : :
2248 : : //=============================================================================
2249 : : //Function: delete_facet (PUBLIC)
2250 : : //Description: delete a single facet and its underlying edges and points if
2251 : : // they are no longer attached to anything
2252 : : //Author: sjowen
2253 : : //Date: 1/21/2004
2254 : : //=============================================================================
2255 : 0 : void FacetDataUtil::delete_facet(CubitFacet *facet_ptr)
2256 : : {
2257 [ # # ]: 0 : DLIList<CubitPoint *>point_list;
2258 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *>edge_list;
2259 [ # # ]: 0 : facet_ptr->points(point_list);
2260 [ # # ]: 0 : facet_ptr->edges(edge_list);
2261 : :
2262 [ # # ][ # # ]: 0 : delete facet_ptr;
2263 : :
2264 : : CubitFacetEdge *edge_ptr;
2265 : : CubitPoint *point_ptr;
2266 : : int ii;
2267 : :
2268 [ # # ][ # # ]: 0 : for (ii=0; ii<edge_list.size(); ii++)
2269 : : {
2270 [ # # ]: 0 : edge_ptr = edge_list.get_and_step();
2271 [ # # ][ # # ]: 0 : if (edge_ptr->num_adj_facets() == 0)
2272 [ # # ][ # # ]: 0 : delete edge_ptr;
2273 : : }
2274 : :
2275 [ # # ]: 0 : for (ii=0; ii<3; ii++)
2276 : : {
2277 [ # # ]: 0 : point_ptr = point_list.get_and_step();
2278 [ # # ][ # # ]: 0 : if (point_ptr->num_adj_facets() == 0)
2279 [ # # ][ # # ]: 0 : delete point_ptr;
2280 [ # # ]: 0 : }
2281 : :
2282 : 0 : }
2283 : :
2284 : :
2285 : 0 : void FacetDataUtil::destruct_facet_no_delete(CubitFacet *facet_ptr)
2286 : : {
2287 [ # # ]: 0 : CubitFacetData* facet_d_ptr = dynamic_cast<CubitFacetData*>(facet_ptr);
2288 [ # # ]: 0 : if(!facet_d_ptr){
2289 [ # # ][ # # ]: 0 : PRINT_ERROR("Can't work with Facet pointer that isn't a facet data object.\n");
[ # # ][ # # ]
2290 : 0 : return;
2291 : : }
2292 : :
2293 [ # # ]: 0 : DLIList<CubitPoint *>point_list;
2294 [ # # ][ # # ]: 0 : DLIList<CubitFacetEdge *>edge_list;
2295 [ # # ]: 0 : facet_ptr->points(point_list);
2296 [ # # ]: 0 : facet_ptr->edges(edge_list);
2297 : :
2298 [ # # ]: 0 : facet_d_ptr->destruct_facet_internals();
2299 : :
2300 : : CubitFacetEdge *edge_ptr;
2301 : : CubitPoint *point_ptr;
2302 : : int ii;
2303 : :
2304 [ # # ][ # # ]: 0 : for (ii=0; ii<edge_list.size(); ii++)
2305 : : {
2306 [ # # ]: 0 : edge_ptr = edge_list.get_and_step();
2307 [ # # ][ # # ]: 0 : if (edge_ptr->num_adj_facets() == 0)
2308 [ # # ][ # # ]: 0 : delete edge_ptr;
2309 : : }
2310 : :
2311 [ # # ]: 0 : for (ii=0; ii<3; ii++)
2312 : : {
2313 [ # # ]: 0 : point_ptr = point_list.get_and_step();
2314 [ # # ][ # # ]: 0 : if (point_ptr->num_adj_facets() == 0)
2315 [ # # ][ # # ]: 0 : delete point_ptr;
2316 [ # # ]: 0 : }
2317 : :
2318 : : }
2319 : :
2320 : : //=============================================================================
2321 : : //Function: intersect_facet (PUBLIC)
2322 : : //Description: determine intersection of a segment with a facet
2323 : : // returns CUBIT_PNT_UNKNOWN: if segment is in plane of facet
2324 : : // CUBIT_PNT_OUTSIDE: if segment does not intersect
2325 : : // CUBIT_PNT_INSIDE: if segment intersects inside facet
2326 : : // CUBIT_PNT_BOUNDARY: if segment intersects a vertex or edge
2327 : : //Author: sjowen
2328 : : //Date: 1/30/2004
2329 : : //=============================================================================
2330 : 0 : CubitPointContainment FacetDataUtil::intersect_facet(
2331 : : CubitVector &start, CubitVector &end, // start and end points of vector
2332 : : CubitFacet *facet_ptr, // the facet to intersect
2333 : : CubitVector &qq, // return the intersection point
2334 : : CubitVector &ac, // area coordinates of qq if is in or on facet
2335 : : CubitBoolean bound) // if true, only check for intersections between the end points.
2336 : : {
2337 : :
2338 [ # # ][ # # ]: 0 : CubitPlane fplane = facet_ptr->plane();
2339 [ # # ]: 0 : double dstart = fplane.distance(start);
2340 [ # # ]: 0 : double dend = fplane.distance(end);
2341 : :
2342 : : // points are both in the plane of the facet - can't handle this case
2343 : :
2344 [ # # ][ # # ]: 0 : if (fabs(dstart) < GEOMETRY_RESABS &&
2345 : 0 : fabs(dend) < GEOMETRY_RESABS)
2346 : : {
2347 : 0 : return CUBIT_PNT_UNKNOWN;
2348 : : }
2349 : :
2350 : : // one point is on the plane
2351 : :
2352 [ # # ]: 0 : if (fabs(dstart) < GEOMETRY_RESABS)
2353 : : {
2354 [ # # ]: 0 : qq = start;
2355 : : }
2356 [ # # ]: 0 : else if (fabs(dend) < GEOMETRY_RESABS)
2357 : : {
2358 [ # # ]: 0 : qq = end;
2359 : : }
2360 : :
2361 : : // points are both on the same side of the plane
2362 [ # # ][ # # ]: 0 : else if(dstart*dend > 0.0 &&
2363 [ # # ]: 0 : (bound || fabs(dstart-dend) < CUBIT_RESABS) )
2364 : : {
2365 : 0 : return CUBIT_PNT_OUTSIDE;
2366 : : }
2367 : : // points are on opposite sides of plane: if bound == false then compute intersection with plane
2368 : :
2369 : : else
2370 : : {
2371 [ # # ]: 0 : CubitVector dir = end-start;
2372 [ # # ]: 0 : dir.normalize();
2373 [ # # ][ # # ]: 0 : qq = fplane.intersect(start, dir);
2374 : : }
2375 : :
2376 [ # # ]: 0 : FacetEvalTool::facet_area_coordinate( facet_ptr, qq, ac );
2377 : :
2378 : : //mod mbrewer ... the original code would call a point
2379 : : // on the boundary if any area coordinate was near
2380 : : // zero, regardless of whether another area coordinate
2381 : : // was negative... making it outside.
2382 : : // if (fabs(ac.x()) < GEOMETRY_RESABS ||
2383 : : // fabs(ac.y()) < GEOMETRY_RESABS ||
2384 : : // fabs(ac.z()) < GEOMETRY_RESABS)
2385 : : // {
2386 : : // return CUBIT_PNT_BOUNDARY;
2387 : : // }
2388 [ # # ][ # # ]: 0 : if ( (fabs(ac.x()) < GEOMETRY_RESABS && (ac.y() > -GEOMETRY_RESABS &&
[ # # ][ # # ]
[ # # ]
2389 [ # # ][ # # ]: 0 : ac.z() > -GEOMETRY_RESABS) )||
2390 [ # # ][ # # ]: 0 : (fabs(ac.y()) < GEOMETRY_RESABS && (ac.x() > -GEOMETRY_RESABS &&
[ # # ][ # # ]
2391 [ # # ][ # # ]: 0 : ac.z() > -GEOMETRY_RESABS) )||
[ # # ]
2392 [ # # ][ # # ]: 0 : (fabs(ac.z()) < GEOMETRY_RESABS && (ac.x() > -GEOMETRY_RESABS &&
[ # # ][ # # ]
2393 [ # # ]: 0 : ac.y() > -GEOMETRY_RESABS) ) ){
2394 : 0 : return CUBIT_PNT_BOUNDARY;
2395 : : }
2396 [ # # ][ # # ]: 0 : else if (ac.x() < 0.0 || ac.y() < 0.0 || ac.z() < 0.0)
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
2397 : : {
2398 : 0 : return CUBIT_PNT_OUTSIDE;
2399 : : }
2400 : :
2401 : 0 : return CUBIT_PNT_INSIDE;
2402 : : }
2403 : :
2404 : :
2405 : : //=============================================================================
2406 : : //Function: get_bbox_of_points (PUBLIC)
2407 : : //Description: Find the bounding box of a list of CubitPoints
2408 : : //Author: jdfowle
2409 : : //Date: 12/15/03
2410 : : //=============================================================================
2411 : 0 : CubitStatus FacetDataUtil::get_bbox_of_points(DLIList<CubitPoint*>& point_list, CubitBox& bbox)
2412 : : {
2413 : : double x, y, z, min[3], max[3];
2414 : : int i;
2415 : : CubitPoint *point;
2416 : 0 : min[0] = min[1] = min[2] = CUBIT_DBL_MAX;
2417 : 0 : max[0] = max[1] = max[2] = -CUBIT_DBL_MAX + 1.;
2418 : :
2419 [ # # ][ # # ]: 0 : for ( i = 0; i < point_list.size(); i++ ) {
2420 [ # # ]: 0 : point = point_list.get_and_step();
2421 [ # # ]: 0 : x = point->x();
2422 [ # # ]: 0 : if ( min[0] > x ) min[0] = x;
2423 [ # # ]: 0 : if ( max[0] < x ) max[0] = x;
2424 [ # # ]: 0 : y = point->y();
2425 [ # # ]: 0 : if ( min[1] > y ) min[1] = y;
2426 [ # # ]: 0 : if ( max[1] < y ) max[1] = y;
2427 [ # # ]: 0 : z = point->z();
2428 [ # # ]: 0 : if ( min[2] > z ) min[2] = z;
2429 [ # # ]: 0 : if ( max[2] < z ) max[2] = z;
2430 : : }
2431 [ # # ]: 0 : bbox.reset(min,max);
2432 : :
2433 : 0 : return CUBIT_SUCCESS;
2434 : : }
2435 : :
2436 : : //=============================================================================
2437 : : //Function: squared_distance_to_segment (PUBLIC)
2438 : : //Description: get square of distance from point to closest point on a line segment;
2439 : : // returns closest point. Taken from "Geometric Tools for Computer Graphics",
2440 : : // by Schneider & Eberly, sec. 6.1.
2441 : : //Author: jdfowle
2442 : : //Date: 03/08/03
2443 : : //=============================================================================
2444 : 3740 : CubitVector FacetDataUtil::squared_distance_to_segment(CubitVector &p, CubitVector &p0,
2445 : : CubitVector &p1, double &distance2)
2446 : : {
2447 [ + - ][ + - ]: 3740 : CubitVector YmPO, D;
2448 : : double t;
2449 : :
2450 [ + - ][ + - ]: 3740 : D = p1 - p0;
2451 [ + - ][ + - ]: 3740 : YmPO = p - p0;
2452 [ + - ][ + - ]: 3740 : t = D.x()*YmPO.x() + D.y()*YmPO.y() + D.z()*YmPO.z();
[ + - ][ + - ]
[ + - ][ + - ]
2453 : :
2454 [ + + ]: 3740 : if ( t < 0.0 ) {
2455 [ + - ][ + - ]: 187 : distance2 = YmPO.x()*YmPO.x() + YmPO.y()*YmPO.y() + YmPO.z()*YmPO.z();
[ + - ][ + - ]
[ + - ][ + - ]
2456 [ + - ]: 187 : return p0;
2457 : : }
2458 : :
2459 : : double DdD;
2460 [ + - ][ + - ]: 3553 : DdD = D.x()*D.x() + D.y()*D.y() + D.z()*D.z();
[ + - ][ + - ]
[ + - ][ + - ]
2461 [ + + ]: 3553 : if ( t >= DdD ) {
2462 [ + - ]: 1551 : CubitVector YmP1;
2463 [ + - ][ + - ]: 1551 : YmP1 = p - p1;
2464 [ + - ][ + - ]: 1551 : distance2 = YmP1.x()*YmP1.x() + YmP1.y()*YmP1.y() + YmP1.z()*YmP1.z();
[ + - ][ + - ]
[ + - ][ + - ]
2465 [ + - ]: 1551 : return p1;
2466 : : }
2467 : :
2468 : : // closest point is interior to segment
2469 [ + - ][ + - ]: 2002 : distance2 = YmPO.x()*YmPO.x() + YmPO.y()*YmPO.y() + YmPO.z()*YmPO.z() - t*t/DdD;
[ + - ][ + - ]
[ + - ][ + - ]
2470 [ + - ][ + - ]: 3740 : return p0 + (t/DdD)*(p1 - p0);
[ + - ]
2471 : : }
2472 : :
2473 : : //=============================================================================
2474 : : //Function: get_bbox_intersections (PUBLIC)
2475 : : //Description: Get the intersection of the line defined by point1 and point2 with
2476 : : //bbox. Returns 0,1 or 2 for the number of intersections. A line
2477 : : //in one of the planes of the box will return 0. Returns -1 for failure.
2478 : : //Author mborden
2479 : : //Date: 04/05/07
2480 : : //=============================================================================
2481 : 0 : int FacetDataUtil::get_bbox_intersections(CubitVector& point1,
2482 : : CubitVector& point2,
2483 : : const CubitBox& bbox,
2484 : : CubitVector& intersection_1,
2485 : : CubitVector& intersection_2)
2486 : : {
2487 : 0 : int debug = 0;
2488 [ # # ]: 0 : if( debug )
2489 : : {
2490 [ # # ]: 0 : GfxDebug::draw_point( point1, CUBIT_RED_INDEX );
2491 [ # # ]: 0 : GfxDebug::draw_point( point2, CUBIT_BLUE_INDEX );
2492 [ # # ]: 0 : GfxDebug::flush();
2493 : : }
2494 : :
2495 : : double coords[6];
2496 [ # # ]: 0 : coords[0] = bbox.min_x();
2497 [ # # ]: 0 : coords[1] = bbox.max_x();
2498 [ # # ]: 0 : coords[2] = bbox.min_y();
2499 [ # # ]: 0 : coords[3] = bbox.max_y();
2500 [ # # ]: 0 : coords[4] = bbox.min_z();
2501 [ # # ]: 0 : coords[5] = bbox.max_z();
2502 : :
2503 [ # # ]: 0 : DLIList<CubitVector*> intersections;
2504 : :
2505 : : int ii;
2506 [ # # ]: 0 : for( ii = 0; ii < 3; ii++ )
2507 : : {
2508 : : //Define four points for each plane.
2509 : : double box_points[4][3];
2510 : :
2511 : : //ii = 0 -> x-planes
2512 : : //ii = 1 -> y_planes
2513 : : //ii = 2 -> z_planes
2514 : :
2515 : : //Only the coordinates for the plane we are in
2516 : : //change. The other two are constant for the
2517 : : //jj loops below.
2518 : 0 : box_points[0][(ii + 1) % 3] = coords[((ii*2)+2) % 6];
2519 : 0 : box_points[1][(ii + 1) % 3] = coords[((ii*2)+3) % 6];
2520 : 0 : box_points[2][(ii + 1) % 3] = coords[((ii*2)+3) % 6];
2521 : 0 : box_points[3][(ii + 1) % 3] = coords[((ii*2)+2) % 6];
2522 : :
2523 : 0 : box_points[0][(ii + 2) % 3] = coords[((ii*2)+4) % 6];
2524 : 0 : box_points[1][(ii + 2) % 3] = coords[((ii*2)+4) % 6];
2525 : 0 : box_points[2][(ii + 2) % 3] = coords[((ii*2)+5) % 6];
2526 : 0 : box_points[3][(ii + 2) % 3] = coords[((ii*2)+5) % 6];
2527 : :
2528 : : int jj;
2529 [ # # ]: 0 : for( jj = 0; jj < 2; jj++ )
2530 : : {
2531 : : CubitPoint* points[4];
2532 : : int kk;
2533 [ # # ]: 0 : for( kk = 0; kk < 4; kk++ )
2534 : : {
2535 : 0 : box_points[kk][ii] = coords[(ii*2)+jj];
2536 [ # # ][ # # ]: 0 : points[kk] = new CubitPointData( box_points[kk][0], box_points[kk][1], box_points[kk][2] );
2537 : : }
2538 : :
2539 : : //Create two facets for this plane to check for intersections.
2540 : : CubitFacet* facets[2];
2541 [ # # ][ # # ]: 0 : facets[0] = new CubitFacetData( points[0], points[1], points[3] );
2542 [ # # ][ # # ]: 0 : facets[1] = new CubitFacetData( points[1], points[2], points[3] );
2543 : :
2544 [ # # ]: 0 : for( kk = 0; kk < 2; kk++ )
2545 : : {
2546 [ # # ]: 0 : CubitVector intersection;
2547 [ # # ]: 0 : CubitVector area_coord;
2548 : :
2549 : : //Make sure the points are not parrellel with the facet.
2550 [ # # ]: 0 : CubitVector dir = point2 - point1;
2551 [ # # ]: 0 : CubitVector normal = facets[kk]->normal();
2552 [ # # ][ # # ]: 0 : if( fabs(dir % normal) < CUBIT_RESABS )
2553 : 0 : continue;
2554 : :
2555 : : CubitPointContainment contain = intersect_facet( point1, point2, facets[kk],
2556 [ # # ]: 0 : intersection, area_coord, CUBIT_FALSE );
2557 [ # # ]: 0 : if( CUBIT_PNT_UNKNOWN == contain )
2558 : : {
2559 : : //The points are in a plane. Return 0.
2560 [ # # ][ # # ]: 0 : delete facets[0];
2561 [ # # ][ # # ]: 0 : delete facets[1];
2562 : : int ll;
2563 [ # # ]: 0 : for( ll = 0; ll < 4; ll++ )
2564 [ # # ][ # # ]: 0 : delete points[ll];
2565 [ # # ][ # # ]: 0 : for( ll = intersections.size(); ll > 0; ll-- )
2566 [ # # ]: 0 : delete intersections.get_and_step();
2567 : :
2568 : 0 : return 0;
2569 : : }
2570 [ # # ][ # # ]: 0 : if( CUBIT_PNT_BOUNDARY == contain ||
2571 : : CUBIT_PNT_INSIDE == contain )
2572 : : {
2573 : : //The point intersects the facet so it's inside the box's surface.
2574 [ # # ][ # # ]: 0 : CubitVector* new_intersection = new CubitVector;
2575 [ # # ]: 0 : *new_intersection = intersection;
2576 [ # # ]: 0 : intersections.append( new_intersection );
2577 : :
2578 [ # # ]: 0 : if( debug )
2579 : : {
2580 [ # # ]: 0 : GfxDebug::draw_point( *new_intersection, CUBIT_CYAN_INDEX );
2581 [ # # ]: 0 : GfxDebug::flush();
2582 : : }
2583 : :
2584 : 0 : break;
2585 : : }
2586 : : }
2587 : :
2588 [ # # ][ # # ]: 0 : delete facets[0];
2589 [ # # ][ # # ]: 0 : delete facets[1];
2590 [ # # ]: 0 : for( kk = 0; kk < 4; kk++ )
2591 [ # # ][ # # ]: 0 : delete points[kk];
2592 : : }
2593 : : }
2594 : :
2595 : : //Check for duplicate intersections.
2596 [ # # ]: 0 : intersections.reset();
2597 [ # # ][ # # ]: 0 : for( ii = 0; ii < intersections.size(); ii++ )
2598 : : {
2599 [ # # ]: 0 : CubitVector* base_vec = intersections.next(ii);
2600 [ # # ]: 0 : if( NULL == base_vec )
2601 : 0 : continue;
2602 : :
2603 : : int jj;
2604 [ # # ][ # # ]: 0 : for( jj = ii+1; jj < intersections.size(); jj++ )
2605 : : {
2606 [ # # ]: 0 : CubitVector* compare_vec = intersections.next(jj);
2607 [ # # ]: 0 : if( NULL != compare_vec )
2608 : : {
2609 [ # # ][ # # ]: 0 : if( base_vec->distance_between_squared( *compare_vec ) < GEOMETRY_RESABS * GEOMETRY_RESABS )
2610 : : {
2611 [ # # ]: 0 : intersections.step(jj);
2612 [ # # ]: 0 : delete intersections.get();
2613 [ # # ]: 0 : intersections.change_to( NULL );
2614 [ # # ]: 0 : intersections.reset();
2615 : : }
2616 : : }
2617 : : }
2618 : : }
2619 [ # # ]: 0 : intersections.remove_all_with_value( NULL );
2620 : :
2621 : :
2622 [ # # ][ # # ]: 0 : if( intersections.size() > 2 )
2623 : : {
2624 [ # # ][ # # ]: 0 : assert( intersections.size() <= 2 );
2625 : 0 : return -1;
2626 : : }
2627 [ # # ][ # # ]: 0 : else if( intersections.size() > 0 )
2628 : : {
2629 [ # # ][ # # ]: 0 : intersection_1 = *intersections.get();
2630 [ # # ][ # # ]: 0 : if( intersections.size() > 1 )
2631 [ # # ][ # # ]: 0 : intersection_2 = *intersections.next();
2632 : : }
2633 : :
2634 : : //Delete memory.
2635 [ # # ][ # # ]: 0 : for( ii = intersections.size(); ii > 0; ii-- )
2636 [ # # ]: 0 : delete intersections.get_and_step();
2637 : :
2638 [ # # ][ # # ]: 0 : return intersections.size();
2639 : : }
2640 : :
2641 : :
2642 : : //=============================================================================
2643 : : //Function: mark_facets (PUBLIC)
2644 : : //Description: mark facets and their children. assumes facets have points and edges
2645 : : //Author sjowen
2646 : : //Date: 09/18/09
2647 : : //=============================================================================
2648 : 0 : void FacetDataUtil::mark_facets( DLIList<FacetEntity *> &facet_list, int mark_value )
2649 : : {
2650 : : int ifacet;
2651 : : FacetEntity *facet_ptr;
2652 [ # # ]: 0 : for (ifacet = 0; ifacet<facet_list.size(); ifacet++ )
2653 : : {
2654 : 0 : facet_ptr = facet_list.get_and_step();
2655 [ # # ]: 0 : CubitFacet *cfacet_ptr = dynamic_cast<CubitFacet *> (facet_ptr);
2656 [ # # ]: 0 : if( cfacet_ptr )
2657 : : {
2658 : 0 : cfacet_ptr->marked(mark_value);
2659 [ # # ]: 0 : for (int ii=0; ii<3; ii++)
2660 : : {
2661 : 0 : cfacet_ptr->point(ii)->marked(mark_value);
2662 : 0 : cfacet_ptr->edge(ii)->marked(mark_value);
2663 : : }
2664 : : }
2665 : : else
2666 : : {
2667 [ # # ]: 0 : CubitFacetEdge *edge_ptr = dynamic_cast<CubitFacetEdge*>(facet_ptr);
2668 [ # # ]: 0 : if( edge_ptr )
2669 : : {
2670 : 0 : edge_ptr->marked(mark_value);
2671 [ # # ]: 0 : for (int ii=0; ii<2; ii++)
2672 : 0 : edge_ptr->point(ii)->marked(mark_value);
2673 : : }
2674 : : }
2675 : : }
2676 : 0 : }
2677 : :
2678 : : //==================================================================================
2679 : : // Description: identifies all the facets that contain given two cubit points
2680 : : //
2681 : : //
2682 : : // Notes:
2683 : : // Author: william roshan quadros
2684 : : // Date: 3/30/2011
2685 : : //===================================================================================
2686 : 0 : CubitStatus FacetDataUtil::find_facet( DLIList<CubitFacet *> temp_split_facets, CubitPoint *pnt0, CubitPoint *pnt1, DLIList<CubitFacet *> &select_facets )
2687 : : {
2688 : : int i;
2689 : : CubitFacet *ptr_facet;
2690 [ # # ][ # # ]: 0 : for( i = 0; i < temp_split_facets.size(); i++ )
2691 : : {
2692 [ # # ]: 0 : ptr_facet = temp_split_facets.get_and_step();
2693 : :
2694 [ # # ][ # # ]: 0 : if( ptr_facet->contains( pnt0 ) && ptr_facet->contains( pnt1 ) )
[ # # ][ # # ]
[ # # ]
2695 : : {
2696 [ # # ]: 0 : select_facets.append( ptr_facet);
2697 : : }
2698 : : }
2699 : :
2700 [ # # ][ # # ]: 0 : if( select_facets.size() )
2701 : 0 : return CUBIT_SUCCESS;
2702 : : else
2703 : 0 : return CUBIT_FAILURE;
2704 [ + - ][ + - ]: 6540 : }
2705 : :
2706 : :
2707 : :
|
