Branch data Line data Source code
1 : : /*
2 : : *
3 : : *
4 : : * Copyright (C) 2004 Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000
5 : : * with Sandia Corporation, the U.S. Government retains certain rights in this software.
6 : : *
7 : : * This file is part of facetbool--contact via [email protected]
8 : : *
9 : : * This library is free software; you can redistribute it and/or
10 : : * modify it under the terms of the GNU Lesser General Public
11 : : * License as published by the Free Software Foundation; either
12 : : * version 2.1 of the License, or (at your option) any later version.
13 : : *
14 : : * This library is distributed in the hope that it will be useful,
15 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : * Lesser General Public License for more details.
18 : : *
19 : : * You should have received a copy of the GNU Lesser General Public
20 : : * License along with this library; if not, write to the Free Software
21 : : * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : : *
23 : : *
24 : : *
25 : : */
26 : :
27 : : #include "FBDataUtil.hpp"
28 : : #include "GeometryDefines.h"
29 : : #include "DLIList.hpp"
30 : :
31 : : //===========================================================================
32 : : //Function Name: intersect_plane_with_boundingbox
33 : : //Description: Find and return the number of intersections of a plane with
34 : : // the edges of a bounding box. There will be 6 or fewer.
35 : : // These will be returned in the intersection_points array.
36 : : //Author: mlstate - rewritten from John Fowler's original version.
37 : : //Date: 3/2005
38 : : //===========================================================================
39 : 0 : void FBDataUtil::intersect_plane_with_boundingbox(CubitBox &bbox,
40 : : const CubitVector &v0,
41 : : const CubitVector &v1,
42 : : const CubitVector &v2,
43 : : DLIList<CubitVector> &intersection_points )
44 : : {
45 [ # # ]: 0 : CubitVector v3 = v0 - v1;
46 [ # # ]: 0 : CubitVector v4 = v2 - v1;
47 [ # # ]: 0 : CubitVector threepointnormal = v3*v4;
48 [ # # ]: 0 : threepointnormal.normalize();
49 : : //CubitVector threepointcenter = (v0 + v1 + v2)/3.;
50 : : // Need to get the size and location of the facetted cutting plane. Do this by
51 : : // intersecting the 3-point plane with the bounding box edges. Will get a max of
52 : : // six intersections. The size of the box for these intersections will determine
53 : : // the cutting plane size; the centroid of the intersections will be the location
54 : : // of the center of the cutting plane.
55 : :
56 : : double a, b, c, d; // plane coefficients
57 [ # # ]: 0 : a = threepointnormal.x();
58 [ # # ]: 0 : b = threepointnormal.y();
59 [ # # ]: 0 : c = threepointnormal.z();
60 [ # # ][ # # ]: 0 : d = -(a*v0.x() + b*v0.y() + c*v0.z());
[ # # ]
61 : : double xbmin, ybmin, zbmin, xbmax, ybmax, zbmax;
62 : : double testpoint;
63 : :
64 [ # # ]: 0 : CubitVector boxmin = bbox.minimum();
65 [ # # ][ # # ]: 0 : xbmin = boxmin.x(); ybmin = boxmin.y(); zbmin = boxmin.z();
[ # # ]
66 [ # # ]: 0 : CubitVector boxmax = bbox.maximum();
67 [ # # ][ # # ]: 0 : xbmax = boxmax.x(); ybmax = boxmax.y(); zbmax = boxmax.z();
[ # # ]
68 : :
69 [ # # ][ # # ]: 0 : if ( fabs(a) < GEOMETRY_RESABS &&
70 : 0 : fabs(b) < GEOMETRY_RESABS )
71 : : {
72 : : // Normal points in Z dir only.
73 : : // a = 0; b = 0
74 [ # # ][ # # ]: 0 : if ( ((zbmin + d/c) < -GEOMETRY_RESABS) && ((zbmax + d/c) > GEOMETRY_RESABS) )
75 : : {
76 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmin, -d/c) );
77 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmin, -d/c) );
78 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmax, -d/c) );
79 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmax, -d/c) );
80 : : }
81 : : }
82 [ # # ][ # # ]: 0 : else if ( fabs(a) < GEOMETRY_RESABS &&
83 : 0 : fabs(c) < GEOMETRY_RESABS )
84 : : {
85 : : // Normal points in Y dir only.
86 [ # # ][ # # ]: 0 : if ( ((ybmin + d/b) < -GEOMETRY_RESABS) && ((ybmax + d/b) > GEOMETRY_RESABS) )
87 : : {
88 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, -d/b, zbmin ) );
89 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, -d/b, zbmin ) );
90 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, -d/b, zbmax ) );
91 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, -d/b, zbmax ) );
92 : : }
93 : : }
94 [ # # ][ # # ]: 0 : else if ( fabs(b) < GEOMETRY_RESABS &&
95 : 0 : fabs(c) < GEOMETRY_RESABS )
96 : : {
97 : : // Normal points in X dir only.
98 [ # # ][ # # ]: 0 : if ( ((xbmin + d/a) < -GEOMETRY_RESABS) && ((xbmax + d/a) > GEOMETRY_RESABS) )
99 : : {
100 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(-d/a, ybmin, zbmin ) );
101 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(-d/a, ybmax, zbmin ) );
102 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(-d/a, ybmax, zbmax ) );
103 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(-d/a, ybmin, zbmax ) );
104 : : }
105 : : }
106 [ # # ]: 0 : else if ( fabs(a) < GEOMETRY_RESABS )
107 : : {
108 : : // Normal is in YZ plane.
109 : 0 : testpoint = -(c*zbmin + d)/b;
110 [ # # ][ # # ]: 0 : if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
111 : 0 : (ybmax > (testpoint - GEOMETRY_RESABS)) )
112 : : {
113 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
114 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
115 : : }
116 : 0 : testpoint = -(c*zbmax + d)/b;
117 [ # # ][ # # ]: 0 : if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
118 : 0 : (ybmax > (testpoint - GEOMETRY_RESABS)) )
119 : : {
120 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
121 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
122 : : }
123 : 0 : testpoint = -(b*ybmin + d)/c;
124 [ # # ][ # # ]: 0 : if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
125 : 0 : (zbmax > (testpoint - GEOMETRY_RESABS)) )
126 : : {
127 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
128 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
129 : : }
130 : 0 : testpoint = -(b*ybmax + d)/c;
131 [ # # ][ # # ]: 0 : if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
132 : 0 : (zbmax > (testpoint - GEOMETRY_RESABS)) )
133 : : {
134 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
135 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
136 : : }
137 : : }
138 [ # # ]: 0 : else if ( fabs(b) < GEOMETRY_RESABS )
139 : : {
140 : : // Normal is in XZ plane
141 : 0 : testpoint = -(c*zbmin + d)/a;
142 [ # # ][ # # ]: 0 : if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
143 : 0 : (xbmax > (testpoint - GEOMETRY_RESABS)) )
144 : : {
145 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
146 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
147 : : }
148 : 0 : testpoint = -(c*zbmax + d)/a;
149 [ # # ][ # # ]: 0 : if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
150 : 0 : (xbmax > (testpoint - GEOMETRY_RESABS)) )
151 : : {
152 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
153 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
154 : : }
155 : 0 : testpoint = -(a*xbmin + d)/c;
156 [ # # ][ # # ]: 0 : if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
157 : 0 : (zbmax > (testpoint - GEOMETRY_RESABS)) )
158 : : {
159 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
160 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
161 : : }
162 : 0 : testpoint = -(a*xbmax + d)/c;
163 [ # # ][ # # ]: 0 : if ( (zbmin < (testpoint + GEOMETRY_RESABS)) &&
164 : 0 : (zbmax > (testpoint - GEOMETRY_RESABS)) )
165 : : {
166 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
167 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
168 : : }
169 : : }
170 [ # # ]: 0 : else if ( fabs(c) < GEOMETRY_RESABS )
171 : : {
172 : : // Normal is in XY plane
173 : 0 : testpoint = -(a*xbmin + d)/b;
174 [ # # ][ # # ]: 0 : if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
175 : 0 : (ybmax > (testpoint - GEOMETRY_RESABS)) )
176 : : {
177 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
178 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
179 : : }
180 : 0 : testpoint = -(a*xbmax + d)/b;
181 [ # # ][ # # ]: 0 : if ( (ybmin < (testpoint + GEOMETRY_RESABS)) &&
182 : 0 : (ybmax > (testpoint - GEOMETRY_RESABS)) )
183 : : {
184 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
185 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
186 : : }
187 : 0 : testpoint = -(b*ybmin + d)/a;
188 [ # # ][ # # ]: 0 : if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
189 : 0 : (xbmax > (testpoint - GEOMETRY_RESABS)) )
190 : : {
191 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
192 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
193 : : }
194 : 0 : testpoint = -(b*ybmax + d)/a;
195 [ # # ][ # # ]: 0 : if ( (xbmin < (testpoint + GEOMETRY_RESABS)) &&
196 : 0 : (xbmax > (testpoint - GEOMETRY_RESABS)) )
197 : : {
198 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
199 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
200 : : }
201 : : }
202 : : else
203 : : {
204 : : // The general case
205 : : // a != 0; b != 0; c != 0
206 : 0 : testpoint = -(b*ybmin + c*zbmin + d)/a;
207 [ # # ][ # # ]: 0 : if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
208 : : {
209 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmin) );
210 : : }
211 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
212 : 0 : testpoint = -(b*ybmax + c*zbmin + d)/a;
213 [ # # ][ # # ]: 0 : if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
214 : : {
215 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmin) );
216 : : }
217 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
218 : 0 : testpoint = -(b*ybmax + c*zbmax + d)/a;
219 [ # # ][ # # ]: 0 : if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
220 : : {
221 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmax, zbmax) );
222 : : }
223 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
224 : 0 : testpoint = -(b*ybmin + c*zbmax + d)/a;
225 [ # # ][ # # ]: 0 : if ( (testpoint > (xbmin-GEOMETRY_RESABS)) && (testpoint < (xbmax+GEOMETRY_RESABS)) )
226 : : {
227 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(testpoint, ybmin, zbmax) );
228 : : }
229 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
230 : 0 : testpoint = -(a*xbmin + c*zbmin + d)/b;
231 [ # # ][ # # ]: 0 : if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
232 : : {
233 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmin) );
234 : : }
235 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
236 : 0 : testpoint = -(a*xbmax + c*zbmin + d)/b;
237 [ # # ][ # # ]: 0 : if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
238 : : {
239 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmin) );
240 : : }
241 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
242 : 0 : testpoint = -(a*xbmax + c*zbmax + d)/b;
243 [ # # ][ # # ]: 0 : if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
244 : : {
245 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, testpoint, zbmax) );
246 : : }
247 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
248 : 0 : testpoint = -(a*xbmin + c*zbmax + d)/b;
249 [ # # ][ # # ]: 0 : if ( (testpoint > (ybmin-GEOMETRY_RESABS)) && (testpoint < (ybmax+GEOMETRY_RESABS)) )
250 : : {
251 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, testpoint, zbmax) );
252 : : }
253 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
254 : 0 : testpoint = -(a*xbmin + b*ybmin + d)/c;
255 [ # # ][ # # ]: 0 : if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
256 : : {
257 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmin, testpoint) );
258 : : }
259 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
260 : 0 : testpoint = -(a*xbmax + b*ybmin + d)/c;
261 [ # # ][ # # ]: 0 : if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
262 : : {
263 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmin, testpoint) );
264 : : }
265 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
266 : 0 : testpoint = -(a*xbmax + b*ybmax + d)/c;
267 [ # # ][ # # ]: 0 : if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
268 : : {
269 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmax, ybmax, testpoint) );
270 : : }
271 [ # # ][ # # ]: 0 : if ( intersection_points.size() == 6 ) return;
272 : 0 : testpoint = -(a*xbmin + b*ybmax + d)/c;
273 [ # # ][ # # ]: 0 : if ( (testpoint > (zbmin-GEOMETRY_RESABS)) && (testpoint < (zbmax+GEOMETRY_RESABS)) )
274 : : {
275 [ # # ][ # # ]: 0 : add_unique_point( intersection_points, CubitVector(xbmin, ybmax, testpoint) );
276 : : }
277 : : }
278 : : }
279 : :
280 : : //===========================================================================
281 : : //Function Name: add_unique_point
282 : : //Description: Given a point and a list of points, add the point to the list
283 : : // unless it is already in the list within GEOMETRY_RESABS.
284 : : //Author:
285 : : //Date: 3/2005
286 : : //===========================================================================
287 : 0 : void FBDataUtil::add_unique_point
288 : : (
289 : : DLIList<CubitVector > &points,
290 : : const CubitVector &pt
291 : : )
292 : : {
293 : : int ipt;
294 [ # # ]: 0 : for ( ipt = 0; ipt < points.size(); ipt++ )
295 : : {
296 : 0 : double dist = pt.distance_between( points[ipt] );
297 [ # # ]: 0 : if ( dist <= GEOMETRY_RESABS )
298 : : {
299 : 0 : return;
300 : : }
301 : : }
302 : 0 : points.append( pt );
303 : : }
304 : :
305 : : //===========================================================================
306 : : //Function Name: FBmake_xy_plane
307 : : //Description: Makes a triangle plane centered on the origin in the x-y
308 : : // direction.
309 : : //Author: jdfowle
310 : : //Date: 1/2004
311 : : //===========================================================================
312 : 0 : CubitStatus FBDataUtil::FBmake_xy_plane(std::vector<double>& verts,
313 : : std::vector<int>& conns,
314 : : double xsize,
315 : : double ysize,
316 : : int numx,
317 : : int numy)
318 : : {
319 : : int i, j;
320 : : double xc, yc, zc, xinc, yinc;
321 : :
322 : 0 : zc = 0.;
323 : 0 : xinc = xsize/(double)(numx-1);
324 : 0 : yinc = ysize/(double)(numy-1);
325 : 0 : xc = -0.5*xsize;
326 : : // Make the coordinates
327 [ # # ]: 0 : for ( i = 0; i < numx; i++ ) {
328 : 0 : yc = -0.5*ysize;
329 [ # # ]: 0 : for ( j = 0; j < numy; j++ ) {
330 [ # # ]: 0 : verts.push_back(xc);
331 [ # # ]: 0 : verts.push_back(yc);
332 [ # # ]: 0 : verts.push_back(zc);
333 : 0 : yc += yinc;
334 : : }
335 : 0 : xc += xinc;
336 : : }
337 : : // Make the connections
338 : : int i1, i2, i3, i4;
339 : :
340 [ # # ]: 0 : for ( i = 0; i < numx-1; i++ ) {
341 : :
342 [ # # ]: 0 : for ( j = 0; j < numy-1; j++ ) {
343 : 0 : i1 = numy*i + j;
344 : 0 : i2 = numy*(i+1) + j;
345 : 0 : i3 = numy*i + j + 1;
346 : 0 : i4 = numy*(i+1) + j + 1;
347 [ # # ]: 0 : conns.push_back(i1);
348 [ # # ]: 0 : conns.push_back(i2);
349 [ # # ]: 0 : conns.push_back(i4);
350 [ # # ]: 0 : conns.push_back(i1);
351 [ # # ]: 0 : conns.push_back(i4);
352 [ # # ]: 0 : conns.push_back(i3);
353 : : }
354 : : }
355 : :
356 : 0 : return CUBIT_SUCCESS;
357 : : }
358 : :
359 : : //===========================================================================
360 : : //Function Name: rotate_FB_object
361 : : //Description: Rotates a facetedboolean object so that the z-direction
362 : : // points in the direction of NormalDir
363 : : //Author: jdfowle
364 : : //Date: 1/2004
365 : : //===========================================================================
366 : 0 : CubitStatus FBDataUtil::rotate_FB_object(std::vector<double>& verts,
367 : : CubitVector &NormalDir, CubitVector &CenterPt)
368 : : {
369 : : unsigned int i;
370 : : double l1, l2, l3, m1, m2, m3, n1, n2, n3;
371 : : double tnx, tny, tnz, temp;
372 : : double tx, ty, tz;
373 : : double xpcen, ypcen, zpcen;
374 : :
375 : 0 : tnx = NormalDir.x();
376 : 0 : tny = NormalDir.y();
377 : 0 : tnz = NormalDir.z();
378 : 0 : xpcen = CenterPt.x();
379 : 0 : ypcen = CenterPt.y();
380 : 0 : zpcen = CenterPt.z();
381 : :
382 : 0 : n1 = tnx; n2 = tny; n3 = tnz;
383 : 0 : l1 = n3; l2 = 0.; l3 = -n1;
384 : 0 : m1 = -n1*n2; m2 = n1*n1 + n3*n3; m3 = -n2*n3;
385 : 0 : temp = sqrt(l1*l1 + l3*l3);
386 [ # # ]: 0 : if ( fabs(temp) > 1.e-6 ) {
387 : 0 : l1 /= temp; l3 /= temp;
388 : 0 : temp = sqrt(m1*m1 + m2*m2 + m3*m3);
389 : 0 : m1 /= temp; m2 /= temp; m3 /= temp;
390 : : } else {
391 : 0 : l1 = 1.; l2 = l3 = 0.;
392 : 0 : m1 = m2 = 0.; m3 = 1.;
393 : 0 : n1 = n3 = 0.; n2 = -1.;
394 : : }
395 : : // Rotate the plane, whiuch is normal to the Z axis and through the origin,
396 : : // so that it will be normal to threepointnormal. Also translate it so
397 : : // that the center point, (0,0,0), moves to threepointcenter.
398 [ # # ]: 0 : for ( i = 0; i < verts.size(); i += 3 ) {
399 : 0 : tx = verts[i];
400 : 0 : ty = verts[i+1];
401 : 0 : tz = verts[i+2];
402 : 0 : verts[i] = tx*l1 + ty*m1 + tz*n1 + xpcen;
403 : 0 : verts[i+1] = tx*l2 + ty*m2 + tz*n2 + ypcen;
404 : 0 : verts[i+2] = tx*l3 + ty*m3 + tz*n3 + zpcen;
405 : : }
406 : :
407 : 0 : return CUBIT_SUCCESS;
408 : : }
409 : :
410 : 0 : int FBDataUtil::makeahashvaluefrom_coord(double x, double y, double z, int numhashbins)
411 : : {
412 : : double sum;
413 : :
414 [ # # ]: 0 : if ( fabs(x) < 1.e-3 ) x = 0.0;
415 [ # # ]: 0 : if ( fabs(y) < 1.e-3 ) y = 0.0;
416 [ # # ]: 0 : if ( fabs(z) < 1.e-3 ) z = 0.0;
417 : 0 : sum = (int)(10000.0*fabs(x) + 0.5) +
418 : 0 : (int)(10000.0*fabs(y) + 0.5) +
419 : 0 : (int)(10000.0*fabs(z) + 0.5);
420 : :
421 : 0 : return (int)(sum) % numhashbins;
422 : : }
423 : :
424 : 0 : CubitStatus FBDataUtil::FBmake_cylinder(std::vector<double>& verts,
425 : : std::vector<int>& coords,
426 : : double radius,
427 : : double length,
428 : : int nr,
429 : : int nz)
430 : : {
431 : : int i, j;
432 : : CubitStatus status;
433 : : double cfac, rinc, linc;
434 : : double x, y, z;
435 : : int istart, iend, V3, pend;
436 : : double zoffset, lpos, rpos, xrad, yrad;
437 : :
438 : 0 : status = CUBIT_SUCCESS;
439 : :
440 : 0 : rinc = 360.0/(double)nr;
441 : 0 : linc = length/(double)nz;
442 : 0 : cfac = CUBIT_PI/180.;
443 : :
444 : 0 : istart = 0; iend = nz+1;
445 : 0 : V3 = (nz+1)*nr;
446 : 0 : pend = nz;
447 : :
448 : :
449 : : // Make the points.
450 : :
451 : 0 : zoffset = 0.0;
452 : 0 : lpos = -0.5*length;
453 : 0 : xrad = radius;
454 : 0 : yrad = radius;
455 [ # # ]: 0 : for ( i = istart; i < iend; i++ ) {
456 : 0 : rpos = 10.0;
457 [ # # ]: 0 : for ( j = 0; j < nr; j++ ) {
458 : 0 : x = xrad*cos(cfac*rpos);
459 : 0 : y = yrad*sin(cfac*rpos);
460 : 0 : z = lpos;
461 [ # # ]: 0 : verts.push_back(x);
462 [ # # ]: 0 : verts.push_back(y);
463 [ # # ]: 0 : verts.push_back(z);
464 : 0 : rpos += rinc;
465 : : }
466 : 0 : lpos += linc;
467 : 0 : zoffset += linc;
468 : : }
469 : : // Add the two points on the axis at the ends.
470 [ # # ]: 0 : verts.push_back(0.);
471 [ # # ]: 0 : verts.push_back(0.);
472 [ # # ]: 0 : verts.push_back(-0.5*length);
473 [ # # ]: 0 : verts.push_back(0.);
474 [ # # ]: 0 : verts.push_back(0.);
475 [ # # ]: 0 : verts.push_back(0.5*length);
476 : :
477 : : // Make the triangles.
478 : : int vertnum;
479 : 0 : vertnum = 0;
480 [ # # ]: 0 : for ( i = 0; i < pend; i++ ) {
481 [ # # ]: 0 : for ( j = 0; j < nr-1; j++ ) {
482 : : // facet_ptr = new CubitFacetData( points[vertnum+j],points[vertnum+j+1], points[vertnum+j+nr] );
483 [ # # ]: 0 : coords.push_back(vertnum+j);
484 [ # # ]: 0 : coords.push_back(vertnum+j+1);
485 [ # # ]: 0 : coords.push_back(vertnum+j+nr);
486 : :
487 : : // facet_ptr = new CubitFacetData( points[vertnum+j+1],points[vertnum+j+1+nr], points[vertnum+j+nr] );
488 [ # # ]: 0 : coords.push_back(vertnum+j+1);
489 [ # # ]: 0 : coords.push_back(vertnum+j+1+nr);
490 [ # # ]: 0 : coords.push_back(vertnum+j+nr);
491 : :
492 : : }
493 : : // facet_ptr = new CubitFacetData( points[vertnum],points[vertnum+nr], points[vertnum+2*nr-1] );
494 [ # # ]: 0 : coords.push_back(vertnum);
495 [ # # ]: 0 : coords.push_back(vertnum+nr);
496 [ # # ]: 0 : coords.push_back(vertnum+2*nr-1);
497 : :
498 : : // facet_ptr = new CubitFacetData( points[vertnum+nr-1],points[vertnum], points[vertnum+2*nr-1] );
499 [ # # ]: 0 : coords.push_back(vertnum+nr-1);
500 [ # # ]: 0 : coords.push_back(vertnum);
501 [ # # ]: 0 : coords.push_back(vertnum+2*nr-1);
502 : :
503 : 0 : vertnum += nr;
504 : : }
505 : :
506 : : // Endcap(s)
507 [ # # ]: 0 : for ( i = 0; i < nr-1; i++ ) { // top cap
508 : : // facet_ptr = new CubitFacetData( points[vertnum+i],points[vertnum+i+1], points[V3+1] );
509 [ # # ]: 0 : coords.push_back(vertnum+i);
510 [ # # ]: 0 : coords.push_back(vertnum+i+1);
511 [ # # ]: 0 : coords.push_back(V3+1);
512 : :
513 : : }
514 : : // facet_ptr = new CubitFacetData( points[nr-1+vertnum],points[vertnum], points[V3+1] );
515 [ # # ]: 0 : coords.push_back(vertnum+nr-1);
516 [ # # ]: 0 : coords.push_back(vertnum);
517 [ # # ]: 0 : coords.push_back(V3+1);
518 : :
519 : :
520 [ # # ]: 0 : for ( i = 0; i < nr-1; i++ ) { // bottom cap
521 : : // facet_ptr = new CubitFacetData( points[i+1],points[i], points[V3] );
522 [ # # ]: 0 : coords.push_back(i+1);
523 [ # # ]: 0 : coords.push_back(i);
524 [ # # ]: 0 : coords.push_back(V3);
525 : :
526 : : }
527 : : // facet_ptr = new CubitFacetData( points[0],points[nr-1], points[V3] );
528 [ # # ]: 0 : coords.push_back(0);
529 [ # # ]: 0 : coords.push_back(nr-1);
530 [ # # ]: 0 : coords.push_back(V3);
531 : :
532 : 0 : return status;
533 : :
534 : : }
535 : :
536 : 0 : double FBDataUtil::project_point_to_plane(double *point, double a, double b,
537 : : double c, double d,
538 : : double *projected_pt)
539 : : {
540 : : double signed_distance;
541 : :
542 : 0 : signed_distance = point[0]*a + point[1]*b + point[2]*c + d;
543 : 0 : projected_pt[0] = point[0] - signed_distance*a;
544 : 0 : projected_pt[1] = point[1] - signed_distance*b;
545 : 0 : projected_pt[2] = point[2] - signed_distance*c;
546 : :
547 : :
548 : 0 : return signed_distance;
549 : :
550 : : }
551 : :
552 : :
553 : : #define EPS_SEG_SEG 1.e-6
554 : : // From Schneider and Eberly,"Geometric Tools for Computer Graphics",
555 : : // Chap. 10.8, segment to segment
556 : : // Note: if the segments are parallel, *s and *t are undefined; *parallel
557 : : // is true; *sunclipped = 0.0, and *tunclipped is the corresponding value
558 : : // for t.
559 : 0 : double FBDataUtil::closest_seg_seg_dist(double *p0, double *d0, double *p1,
560 : : double *d1, double *s, double *t,
561 : : double *sunclipped, double *tunclipped,
562 : : bool *parallel)
563 : : {
564 : : double ux, uy, uz;
565 : : double a, b, c, d, e, det;
566 : : double snum, sdenom, tnum, tdenom;
567 : :
568 : 0 : *parallel = false;
569 : :
570 : 0 : ux = p0[0] - p1[0];
571 : 0 : uy = p0[1] - p1[1];
572 : 0 : uz = p0[2] - p1[2];
573 : :
574 : 0 : a = d0[0]*d0[0] + d0[1]*d0[1] + d0[2]*d0[2];
575 : 0 : b = d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2];
576 : 0 : c = d1[0]*d1[0] + d1[1]*d1[1] + d1[2]*d1[2];
577 : 0 : d = d0[0]*ux + d0[1]*uy + d0[2]*uz;
578 : 0 : e = d1[0]*ux + d1[1]*uy + d1[2]*uz;
579 : 0 : det = a*c - b*b;
580 : :
581 [ # # ]: 0 : if ( det < EPS_SEG_SEG ) {
582 : 0 : snum = 0.;
583 : 0 : tnum = e;
584 : 0 : tdenom = c;
585 : 0 : sdenom = 1.;
586 : 0 : *sunclipped = snum/sdenom;
587 : 0 : *tunclipped = tnum/tdenom;
588 : 0 : double f = ux*ux + uy*uy + uz*uz;
589 : 0 : *parallel = true;
590 : 0 : return sqrt(f - e*e/c);
591 : : } else {
592 : 0 : snum = b*e - c*d;
593 : 0 : tnum = a*e - b*d;
594 : 0 : sdenom = det;
595 : 0 : *sunclipped = snum/det;
596 : 0 : *tunclipped = tnum/det;
597 : : }
598 : :
599 : :
600 [ # # ]: 0 : if ( snum < 0.0 ) {
601 : 0 : snum = 0.0;
602 : 0 : tnum = e;
603 : 0 : tdenom = c;
604 [ # # ]: 0 : } else if ( snum > det ) {
605 : 0 : snum = det;
606 : 0 : tnum = e + b;
607 : 0 : tdenom = c;
608 : : } else {
609 : 0 : tdenom = det;
610 : : }
611 : :
612 [ # # ]: 0 : if ( tnum < 0.0 ) {
613 : 0 : tnum = 0.0;
614 [ # # ]: 0 : if ( -d < 0.0 ) {
615 : 0 : snum = 0.0;
616 [ # # ]: 0 : } else if ( -d > a ) {
617 : 0 : snum = sdenom;
618 : : } else {
619 : 0 : snum = -d;
620 : 0 : sdenom = a;
621 : : }
622 [ # # ]: 0 : } else if ( tnum > tdenom ) {
623 : 0 : tnum = tdenom;
624 [ # # ]: 0 : if ( (-d + b) < 0.0 ) {
625 : 0 : snum = 0.0;
626 [ # # ]: 0 : } else if ( (-d + b) > a ) {
627 : 0 : snum = sdenom;
628 : : } else {
629 : 0 : snum = -d + b;
630 : 0 : sdenom = a;
631 : : }
632 : : }
633 : :
634 : 0 : *s = snum/sdenom;
635 : 0 : *t = tnum/tdenom;
636 : :
637 : : double vx, vy, vz;
638 : 0 : vx = p0[0] + (*s*d0[0]) - (p1[0] + (*t*d1[0]));
639 : 0 : vy = p0[1] + (*s*d0[1]) - (p1[1] + (*t*d1[1]));
640 : 0 : vz = p0[2] + (*s*d0[2]) - (p1[2] + (*t*d1[2]));
641 : :
642 : 0 : return sqrt(vx*vx + vy*vy + vz*vz);
643 : :
644 : : }
|
