1 /* ----------------------------------------------------------------------
2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3 http://sparta.sandia.gov
4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5 Sandia National Laboratories
6
7 Copyright (2014) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14
15 #include "geometry.h"
16 #include "math_extra.h"
17
18 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
19 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
20
21 #define EPSSQ 1.0e-16
22 #define EPSSQNEG -1.0e-16
23 #define EPSSELF 1.0e-6
24
25 enum{OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN}; // same as Update
26
27 namespace Geometry {
28
29 /* ----------------------------------------------------------------------
30 compute whether line intersects an orthogonal 2d quad cell
31 intersection is defined as
32 any line pt (interior, vertex) in common with
33 any quad pt (interior, edge, vertex)
34 v0,v1 and norm = 2 vertices of line and unit normal vec
35 lo,hi = opposite corner pts of quad
36 return 1 if intersection, else 0
37 ------------------------------------------------------------------------- */
38
line_quad_intersect(double * v0,double * v1,double * norm,double * lo,double * hi)39 int line_quad_intersect(double *v0, double *v1, double *norm,
40 double *lo, double *hi)
41 {
42 int sum,side;
43 double xlo,xhi,ylo,yhi,param;
44 double b[3],e[3],point[3];
45
46 xlo = lo[0];
47 xhi = hi[0];
48 ylo = lo[1];
49 yhi = hi[1];
50
51 // if either of line vertices are inside quad, intersection
52 // use <= and >= so touching quad surface is same as inside it
53 // important to do this test first, b/c whichside test can be epsilon off
54
55 if (v0[0] >= xlo && v0[0] <= xhi && v0[1] >= ylo && v0[1] <= yhi) return 1;
56 if (v1[0] >= xlo && v1[0] <= xhi && v1[1] >= ylo && v1[1] <= yhi) return 1;
57
58 // if all 4 quad pts are on same side of line, no intersection
59
60 sum = whichside(v0,norm,xlo,ylo,0.0);
61 sum += whichside(v0,norm,xhi,ylo,0.0);
62 sum += whichside(v0,norm,xlo,yhi,0.0);
63 sum += whichside(v0,norm,xhi,yhi,0.0);
64
65 if (sum == 4 || sum == -4) return 0;
66
67 // test 4 quad edges for intersection with line
68 // b,e = begin/end of quad edge line segment
69
70 b[0] = xlo; b[1] = ylo; b[2] = 0.0;
71 e[0] = xhi; e[1] = ylo; e[2] = 0.0;
72 if (line_line_intersect(b,e,v0,v1,norm,point,param,side)) return 1;
73
74 b[0] = xhi; b[1] = ylo; b[2] = 0.0;
75 e[0] = xhi; e[1] = yhi; e[2] = 0.0;
76 if (line_line_intersect(b,e,v0,v1,norm,point,param,side)) return 1;
77
78 b[0] = xhi; b[1] = yhi; b[2] = 0.0;
79 e[0] = xlo; e[1] = yhi; e[2] = 0.0;
80 if (line_line_intersect(b,e,v0,v1,norm,point,param,side)) return 1;
81
82 b[0] = xlo; b[1] = yhi; b[2] = 0.0;
83 e[0] = xlo; e[1] = ylo; e[2] = 0.0;
84 if (line_line_intersect(b,e,v0,v1,norm,point,param,side)) return 1;
85
86 return 0;
87 }
88
89 /* ----------------------------------------------------------------------
90 compute any intersection of edges of orthogonal 2d quad cell with a line
91 line interior to quad cell has no intersection
92 v0,v1 and norm = 2 vertices of line and unit normal vec
93 lo,hi = opposite corner pts of quad
94 return 1 if intersection, else 0
95 return xc = intersection point if there is one
96 ------------------------------------------------------------------------- */
97
quad_line_intersect_point(double * v0,double * v1,double * norm,double * lo,double * hi,double * xc)98 int quad_line_intersect_point(double *v0, double *v1, double *norm,
99 double *lo, double *hi, double *xc)
100 {
101 int side;
102 double xlo,xhi,ylo,yhi,param;
103 double b[3],e[3];
104
105 xlo = lo[0];
106 xhi = hi[0];
107 ylo = lo[1];
108 yhi = hi[1];
109
110 // test 4 quad edges for intersection with line
111 // b,e = begin/end of quad edge line segment
112
113 b[0] = xlo; b[1] = ylo; b[2] = 0.0;
114 e[0] = xhi; e[1] = ylo; e[2] = 0.0;
115 if (line_line_intersect(b,e,v0,v1,norm,xc,param,side)) return 1;
116
117 b[0] = xhi; b[1] = ylo; b[2] = 0.0;
118 e[0] = xhi; e[1] = yhi; e[2] = 0.0;
119 if (line_line_intersect(b,e,v0,v1,norm,xc,param,side)) return 1;
120
121 b[0] = xhi; b[1] = yhi; b[2] = 0.0;
122 e[0] = xlo; e[1] = yhi; e[2] = 0.0;
123 if (line_line_intersect(b,e,v0,v1,norm,xc,param,side)) return 1;
124
125 b[0] = xlo; b[1] = yhi; b[2] = 0.0;
126 e[0] = xlo; e[1] = ylo; e[2] = 0.0;
127 if (line_line_intersect(b,e,v0,v1,norm,xc,param,side)) return 1;
128
129 return 0;
130 }
131
132 /* ----------------------------------------------------------------------
133 compute whether line touches iface of orthogonal 2d quad cell
134 touch is defined as
135 line end pt = any face pt (edge, vertex)
136 v0,v1 = 2 vertices of line
137 iface = 0 to 3 = XLO,XHI,YLO,YHI
138 lo,hi = opposite corner pts of quad
139 return 1 if touches, else 0
140 ------------------------------------------------------------------------- */
141
line_touch_quad_face(double * v0,double * v1,int iface,double * lo,double * hi)142 int line_touch_quad_face(double *v0, double *v1, int iface,
143 double *lo, double *hi)
144 {
145 // value = position of face
146
147 int dim = iface / 2;
148 int other = dim ? 0 : 1;
149 double value = iface % 2 ? hi[dim] : lo[dim];
150
151 // check if either line vertex is within face
152
153 if (v0[dim] == value) {
154 if (v0[other] >= lo[other] && v0[other] <= hi[other]) return 1;
155 }
156 if (v1[dim] == value) {
157 if (v1[other] >= lo[other] && v1[other] <= hi[other]) return 1;
158 }
159
160 return 0;
161 }
162
163 /* ----------------------------------------------------------------------
164 compute whether triangle intersects an orthogonal 3d hex cell
165 intersection is defined as
166 any triangle pt (interior, edge, vertex) in common with
167 any hex pt (interior, face, edge, vertex)
168 v0,v1,v2 and norm = 3 vertices of triangle and unit normal vec
169 lo,hi = opposite corner pts of hex
170 return 1 if intersection, else 0
171 ------------------------------------------------------------------------- */
172
tri_hex_intersect(double * v0,double * v1,double * v2,double * norm,double * lo,double * hi)173 int tri_hex_intersect(double *v0, double *v1, double *v2, double *norm,
174 double *lo, double *hi)
175 {
176 int sum,side;
177 double xlo,xhi,ylo,yhi,zlo,zhi,param;
178 double b[3],e[3],h0[3],h1[3],h2[3],h3[3],n[3],point[3];
179
180 xlo = lo[0];
181 xhi = hi[0];
182 ylo = lo[1];
183 yhi = hi[1];
184 zlo = lo[2];
185 zhi = hi[2];
186
187 // if any of 3 tri vertices are inside hex, intersection
188 // use <= and >= so touching hex surface is same as inside it
189 // important to do this test first, b/c whichside test can be epsilon off
190
191 if (v0[0] >= xlo && v0[0] <= xhi && v0[1] >= ylo && v0[1] <= yhi &&
192 v0[2] >= zlo && v0[2] <= zhi) return 1;
193
194 if (v1[0] >= xlo && v1[0] <= xhi && v1[1] >= ylo && v1[1] <= yhi &&
195 v1[2] >= zlo && v1[2] <= zhi) return 1;
196
197 if (v2[0] >= xlo && v2[0] <= xhi && v2[1] >= ylo && v2[1] <= yhi &&
198 v2[2] >= zlo && v2[2] <= zhi) return 1;
199
200 // if all 8 hex pts are on same side of tri plane, no intersection
201
202 sum = whichside(v0,norm,xlo,ylo,zlo);
203 sum += whichside(v0,norm,xhi,ylo,zlo);
204 sum += whichside(v0,norm,xlo,yhi,zlo);
205 sum += whichside(v0,norm,xhi,yhi,zlo);
206 sum += whichside(v0,norm,xlo,ylo,zhi);
207 sum += whichside(v0,norm,xhi,ylo,zhi);
208 sum += whichside(v0,norm,xlo,yhi,zhi);
209 sum += whichside(v0,norm,xhi,yhi,zhi);
210 if (sum == 8 || sum == -8) return 0;
211
212 // test 12 hex edges for intersection with tri
213 // b,e = begin/end of hex edge line segment
214
215 b[0] = xlo; b[1] = ylo; b[2] = zlo;
216 e[0] = xhi; e[1] = ylo; e[2] = zlo;
217 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
218
219 b[0] = xlo; b[1] = yhi; b[2] = zlo;
220 e[0] = xhi; e[1] = yhi; e[2] = zlo;
221 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
222
223 b[0] = xlo; b[1] = ylo; b[2] = zhi;
224 e[0] = xhi; e[1] = ylo; e[2] = zhi;
225 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
226
227 b[0] = xlo; b[1] = yhi; b[2] = zhi;
228 e[0] = xhi; e[1] = yhi; e[2] = zhi;
229 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
230
231 b[0] = xlo; b[1] = ylo; b[2] = zlo;
232 e[0] = xlo; e[1] = yhi; e[2] = zlo;
233 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
234
235 b[0] = xhi; b[1] = ylo; b[2] = zlo;
236 e[0] = xhi; e[1] = yhi; e[2] = zlo;
237 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
238
239 b[0] = xlo; b[1] = ylo; b[2] = zhi;
240 e[0] = xlo; e[1] = yhi; e[2] = zhi;
241 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
242
243 b[0] = xhi; b[1] = ylo; b[2] = zhi;
244 e[0] = xhi; e[1] = yhi; e[2] = zhi;
245 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
246
247 b[0] = xlo; b[1] = ylo; b[2] = zlo;
248 e[0] = xlo; e[1] = ylo; e[2] = zhi;
249 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
250
251 b[0] = xhi; b[1] = ylo; b[2] = zlo;
252 e[0] = xhi; e[1] = ylo; e[2] = zhi;
253 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
254
255 b[0] = xlo; b[1] = yhi; b[2] = zlo;
256 e[0] = xlo; e[1] = yhi; e[2] = zhi;
257 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
258
259 b[0] = xhi; b[1] = yhi; b[2] = zlo;
260 e[0] = xhi; e[1] = yhi; e[2] = zhi;
261 if (line_tri_intersect(b,e,v0,v1,v2,norm,point,param,side)) return 1;
262
263 // test 3 tri edges for intersection with 6 faces of hex
264 // h0,h1,h2,h3 = 4 corner pts of hex face
265 // n = normal to xyz faces, depends on vertex ordering
266 // each face is treated as 2 triangles -> 6 tests per face
267
268 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
269 h1[0] = xlo; h1[1] = yhi; h1[2] = zlo;
270 h2[0] = xlo; h2[1] = yhi; h2[2] = zhi;
271 h3[0] = xlo; h3[1] = ylo; h3[2] = zhi;
272 n[0] = 1.0; n[1] = 0.0; n[2] = 0.0;
273
274 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
275 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
276 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
277 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
278 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
279 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
280
281 h0[0] = h1[0] = h2[0] = h3[0] = xhi;
282
283 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
284 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
285 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
286 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
287 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
288 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
289
290 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
291 h1[0] = xhi; h1[1] = ylo; h1[2] = zlo;
292 h2[0] = xhi; h2[1] = ylo; h2[2] = zhi;
293 h3[0] = xlo; h3[1] = ylo; h3[2] = zhi;
294 n[0] = 0.0; n[1] = -1.0; n[2] = 0.0;
295
296 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
297 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
298 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
299 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
300 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
301 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
302
303 h0[1] = h1[1] = h2[1] = h3[1] = yhi;
304
305 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
306 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
307 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
308 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
309 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
310 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
311
312 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
313 h1[0] = xhi; h1[1] = ylo; h1[2] = zlo;
314 h2[0] = xhi; h2[1] = yhi; h2[2] = zlo;
315 h3[0] = xlo; h3[1] = yhi; h3[2] = zlo;
316 n[0] = 0.0; n[1] = 0.0; n[2] = 1.0;
317
318 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
319 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
320 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
321 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
322 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
323 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
324
325 h0[2] = h1[2] = h2[2] = h3[2] = zhi;
326
327 if (line_tri_intersect(v0,v1,h0,h1,h2,n,point,param,side) ||
328 line_tri_intersect(v1,v2,h0,h1,h2,n,point,param,side) ||
329 line_tri_intersect(v2,v0,h0,h1,h2,n,point,param,side) ||
330 line_tri_intersect(v0,v1,h0,h2,h3,n,point,param,side) ||
331 line_tri_intersect(v1,v2,h0,h2,h3,n,point,param,side) ||
332 line_tri_intersect(v2,v0,h0,h2,h3,n,point,param,side)) return 1;
333
334 return 0;
335 }
336
337 /* ----------------------------------------------------------------------
338 compute any intersection of edges/faces of orthogonal 3d hex cell with a tri
339 tri interior to quad cell has no intersection
340 v0,v1,v2 and norm = 3 vertices of triangle and unit normal vec
341 lo,hi = opposite corner pts of hex
342 return 1 if intersection, else 0
343 return xc = intersection point if there is one
344 ------------------------------------------------------------------------- */
345
hex_tri_intersect_point(double * v0,double * v1,double * v2,double * norm,double * lo,double * hi,double * xc)346 int hex_tri_intersect_point(double *v0, double *v1, double *v2, double *norm,
347 double *lo, double *hi, double *xc)
348 {
349 int side;
350 double xlo,xhi,ylo,yhi,zlo,zhi,param;
351 double b[3],e[3],h0[3],h1[3],h2[3],h3[3],n[3];
352
353 xlo = lo[0];
354 xhi = hi[0];
355 ylo = lo[1];
356 yhi = hi[1];
357 zlo = lo[2];
358 zhi = hi[2];
359
360 // test 12 hex edges for intersection with tri
361 // b,e = begin/end of hex edge line segment
362
363 b[0] = xlo; b[1] = ylo; b[2] = zlo;
364 e[0] = xhi; e[1] = ylo; e[2] = zlo;
365 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
366
367 b[0] = xlo; b[1] = yhi; b[2] = zlo;
368 e[0] = xhi; e[1] = yhi; e[2] = zlo;
369 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
370
371 b[0] = xlo; b[1] = ylo; b[2] = zhi;
372 e[0] = xhi; e[1] = ylo; e[2] = zhi;
373 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
374
375 b[0] = xlo; b[1] = yhi; b[2] = zhi;
376 e[0] = xhi; e[1] = yhi; e[2] = zhi;
377 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
378
379 b[0] = xlo; b[1] = ylo; b[2] = zlo;
380 e[0] = xlo; e[1] = yhi; e[2] = zlo;
381 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
382
383 b[0] = xhi; b[1] = ylo; b[2] = zlo;
384 e[0] = xhi; e[1] = yhi; e[2] = zlo;
385 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
386
387 b[0] = xlo; b[1] = ylo; b[2] = zhi;
388 e[0] = xlo; e[1] = yhi; e[2] = zhi;
389 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
390
391 b[0] = xhi; b[1] = ylo; b[2] = zhi;
392 e[0] = xhi; e[1] = yhi; e[2] = zhi;
393 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
394
395 b[0] = xlo; b[1] = ylo; b[2] = zlo;
396 e[0] = xlo; e[1] = ylo; e[2] = zhi;
397 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
398
399 b[0] = xhi; b[1] = ylo; b[2] = zlo;
400 e[0] = xhi; e[1] = ylo; e[2] = zhi;
401 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
402
403 b[0] = xlo; b[1] = yhi; b[2] = zlo;
404 e[0] = xlo; e[1] = yhi; e[2] = zhi;
405 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
406
407 b[0] = xhi; b[1] = yhi; b[2] = zlo;
408 e[0] = xhi; e[1] = yhi; e[2] = zhi;
409 if (line_tri_intersect(b,e,v0,v1,v2,norm,xc,param,side)) return 1;
410
411 // test 3 tri edges for intersection with 6 faces of hex
412 // h0,h1,h2,h3 = 4 corner pts of hex face
413 // n = normal to xyz faces, depends on vertex ordering
414 // each face is treated as 2 triangles -> 6 tests per face
415
416 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
417 h1[0] = xlo; h1[1] = yhi; h1[2] = zlo;
418 h2[0] = xlo; h2[1] = yhi; h2[2] = zhi;
419 h3[0] = xlo; h3[1] = ylo; h3[2] = zhi;
420 n[0] = 1.0; n[1] = 0.0; n[2] = 0.0;
421
422 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
423 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
424 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
425 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
426 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
427 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
428
429 h0[0] = h1[0] = h2[0] = h3[0] = xhi;
430
431 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
432 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
433 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
434 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
435 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
436 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
437
438 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
439 h1[0] = xhi; h1[1] = ylo; h1[2] = zlo;
440 h2[0] = xhi; h2[1] = ylo; h2[2] = zhi;
441 h3[0] = xlo; h3[1] = ylo; h3[2] = zhi;
442 n[0] = 0.0; n[1] = -1.0; n[2] = 0.0;
443
444 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
445 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
446 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
447 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
448 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
449 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
450
451 h0[1] = h1[1] = h2[1] = h3[1] = yhi;
452
453 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
454 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
455 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
456 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
457 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
458 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
459
460 h0[0] = xlo; h0[1] = ylo; h0[2] = zlo;
461 h1[0] = xhi; h1[1] = ylo; h1[2] = zlo;
462 h2[0] = xhi; h2[1] = yhi; h2[2] = zlo;
463 h3[0] = xlo; h3[1] = yhi; h3[2] = zlo;
464 n[0] = 0.0; n[1] = 0.0; n[2] = 1.0;
465
466 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
467 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
468 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
469 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
470 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
471 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
472
473 h0[2] = h1[2] = h2[2] = h3[2] = zhi;
474
475 if (line_tri_intersect(v0,v1,h0,h1,h2,n,xc,param,side) ||
476 line_tri_intersect(v1,v2,h0,h1,h2,n,xc,param,side) ||
477 line_tri_intersect(v2,v0,h0,h1,h2,n,xc,param,side) ||
478 line_tri_intersect(v0,v1,h0,h2,h3,n,xc,param,side) ||
479 line_tri_intersect(v1,v2,h0,h2,h3,n,xc,param,side) ||
480 line_tri_intersect(v2,v0,h0,h2,h3,n,xc,param,side)) return 1;
481
482 return 0;
483 }
484
485 /* ----------------------------------------------------------------------
486 compute whether triangle touches iface of orthogonal 3d hex cell
487 touch is defined as
488 triangle corner pt = any face pt (interior, edge, vertex)
489 v0,v1,v2 = 3 vertices of triangle
490 iface = 0 to 5 = XLO,XHI,YLO,YHI,ZLO,ZHI
491 lo,hi = opposite corner pts of quad
492 return 1 if touches, else 0
493 ------------------------------------------------------------------------- */
494
tri_touch_hex_face(double * v0,double * v1,double * v2,int iface,double * lo,double * hi)495 int tri_touch_hex_face(double *v0, double *v1, double *v2, int iface,
496 double *lo, double *hi)
497 {
498 // value = position of face
499
500 int dim = iface / 2;
501 int other1,other2;
502 if (dim == 0) {
503 other1 = 1; other2 = 2;
504 } else if (dim == 1) {
505 other1 = 0; other2 = 2;
506 } else if (dim == 2) {
507 other1 = 0; other2 = 1;
508 }
509 double value = iface % 2 ? hi[dim] : lo[dim];
510
511 // check if any triangle vertex is within face
512
513 if (v0[dim] == value) {
514 if (v0[other1] >= lo[other1] && v0[other1] <= hi[other1] &&
515 v0[other2] >= lo[other2] && v0[other2] <= hi[other2]) return 1;
516 }
517 if (v1[dim] == value) {
518 if (v1[other1] >= lo[other1] && v1[other1] <= hi[other1] &&
519 v1[other2] >= lo[other2] && v1[other2] <= hi[other2]) return 1;
520 }
521 if (v2[dim] == value) {
522 if (v2[other1] >= lo[other1] && v2[other1] <= hi[other1] &&
523 v2[other2] >= lo[other2] && v2[other2] <= hi[other2]) return 1;
524 }
525
526 return 0;
527 }
528
529 /* ----------------------------------------------------------------------
530 compute whether triangle sits on any face of orthogonal 3d hex cell
531 return 0 to 5 = XLO,XHI,YLO,YHI,ZLO,ZHI face
532 else return -1
533 ------------------------------------------------------------------------- */
534
tri_on_hex_face(double * v0,double * v1,double * v2,double * lo,double * hi)535 int tri_on_hex_face(double *v0, double *v1, double *v2, double *lo, double *hi)
536 {
537 double vmin[3],vmax[3];
538
539 vmin[0] = MIN(v0[0],v1[0]);
540 vmin[0] = MIN(v2[0],vmin[0]);
541 vmin[1] = MIN(v0[1],v1[1]);
542 vmin[1] = MIN(v2[1],vmin[1]);
543 vmin[2] = MIN(v0[2],v1[2]);
544 vmin[2] = MIN(v2[2],vmin[2]);
545
546 vmax[0] = MAX(v0[0],v1[0]);
547 vmax[0] = MAX(v2[0],vmax[0]);
548 vmax[1] = MAX(v0[1],v1[1]);
549 vmax[1] = MAX(v2[1],vmax[1]);
550 vmax[2] = MAX(v0[2],v1[2]);
551 vmax[2] = MAX(v2[2],vmax[2]);
552
553 if (vmin[0] == vmax[0]) {
554 if (vmin[0] == lo[0]) return 0;
555 if (vmin[0] == hi[0]) return 1;
556 }
557
558 if (vmin[1] == vmax[1]) {
559 if (vmin[1] == lo[1]) return 2;
560 if (vmin[1] == hi[1]) return 3;
561 }
562
563 if (vmin[2] == vmax[2]) {
564 if (vmin[2] == lo[2]) return 4;
565 if (vmin[2] == hi[2]) return 5;
566 }
567
568 return -1;
569 }
570
571 /* ----------------------------------------------------------------------
572 compute whether triangle edge sits on any face of orthogonal 3d hex cell
573 return 0 to 5 = XLO,XHI,YLO,YHI,ZLO,ZHI face
574 else return -1
575 ------------------------------------------------------------------------- */
576
edge_on_hex_face(double * v0,double * v1,double * lo,double * hi)577 int edge_on_hex_face(double *v0, double *v1, double *lo, double *hi)
578 {
579 double vmin[3],vmax[3];
580
581 vmin[0] = MIN(v0[0],v1[0]);
582 vmin[1] = MIN(v0[1],v1[1]);
583 vmin[2] = MIN(v0[2],v1[2]);
584
585 vmax[0] = MAX(v0[0],v1[0]);
586 vmax[1] = MAX(v0[1],v1[1]);
587 vmax[2] = MAX(v0[2],v1[2]);
588
589 if (vmin[0] == vmax[0]) {
590 if (vmin[0] == lo[0]) return 0;
591 if (vmin[0] == hi[0]) return 1;
592 }
593
594 if (vmin[1] == vmax[1]) {
595 if (vmin[1] == lo[1]) return 2;
596 if (vmin[1] == hi[1]) return 3;
597 }
598
599 if (vmin[2] == vmax[2]) {
600 if (vmin[2] == lo[2]) return 4;
601 if (vmin[2] == hi[2]) return 5;
602 }
603
604 return -1;
605 }
606
607 /* ----------------------------------------------------------------------
608 detect intersection between a directed line segment A and line segment B
609 intersection is defined as any A pt (including end pts)
610 in common with any B pt (interior,vertex)
611 one exception is if both A end pts are on infinite line B,
612 then is NOT an intersection
613 start,stop = end points of directed line segment A
614 A must have non-zero length
615 v0,v1 = 2 vertices of line segment B
616 norm = unit vector normal to line segment B
617 pointing OUTSIDE via right-hand rule
618 return TRUE if there is an intersection, else FALSE
619 if TRUE also return:
620 point = pt of intersection
621 param = intersection pt is this fraction along line A (0-1 inclusive)
622 side = side of B that was hit = OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN
623 ------------------------------------------------------------------------- */
624
line_line_intersect(double * start,double * stop,double * v0,double * v1,double * norm,double * point,double & param,int & side,int)625 bool line_line_intersect(double *start, double *stop,
626 double *v0, double *v1, double *norm,
627 double *point, double ¶m, int &side, int)
628 {
629 double vec[3],start2stop[3],edge[3],pvec[3];
630
631 // if start,stop are on same side of line B, no intersection
632 // if start,stop are both on infinite line B, no intersection
633
634 MathExtra::sub3(start,v0,vec);
635 double dotstart = MathExtra::dot3(norm,vec);
636 MathExtra::sub3(stop,v0,vec);
637 double dotstop = MathExtra::dot3(norm,vec);
638
639 if (dotstart < 0.0 && dotstop < 0.0) return false;
640 if (dotstart > 0.0 && dotstop > 0.0) return false;
641 if (dotstart == 0.0 && dotstop == 0.0) return false;
642
643 // param = parametric distance from start to stop
644 // at which line B is intersected
645 // param = must be 0.0 to 1.0 inclusive, else no intersection
646
647 MathExtra::sub3(v0,start,vec);
648 MathExtra::sub3(stop,start,start2stop);
649 param = MathExtra::dot3(norm,vec) / MathExtra::dot3(norm,start2stop);
650
651 if (param < 0.0 || param > 1.0) return false;
652
653 // point = intersection pt with line B
654
655 point[0] = start[0] + param * start2stop[0];
656 point[1] = start[1] + param * start2stop[1];
657 point[2] = 0.0;
658
659 // test if intersection pt is inside line B
660 // edge = line B vector from v0 to v1
661 // pvec = vector from either line B vertex to intersection point
662 // if dot product of edge with pvec < or > 0.0 for each pvec,
663 // intersection point is outside line B
664 // use EPSSQ and EPSSQNEG instead of 0.0 for following case:
665 // intersection pt is on line end pt where 2 lines come together
666 // want it to detect collision with at least one of lines
667 // point can be epsilon away from end pt
668 // this leads to pvec being epsilon vec in opposite dirs for 2 lines
669 // this can lead to dot3() being negative espilon^2 for both lines,
670 // depending on direction of 2 lines
671 // thus this can lead to no collision with either line
672 // typical observed dot values were 1.0e-18, so use EPSSQ = 1.0e-16
673
674 MathExtra::sub3(v1,v0,edge);
675 MathExtra::sub3(point,v0,pvec);
676 if (MathExtra::dot3(edge,pvec) < EPSSQNEG) return false;
677 MathExtra::sub3(point,v1,pvec);
678 if (MathExtra::dot3(edge,pvec) > EPSSQ) return false;
679
680 // there is a valid intersection with line B
681 // set side to ONSURF, OUTSIDE, or INSIDE
682 // if start point is inside or outside then side = same
683 // if particle started on line B, side = ONSURF OUT/IN based on dotstop
684
685 if (dotstart < 0.0) side = INSIDE;
686 else if (dotstart > 0.0) side = OUTSIDE;
687 else if (dotstop > 0.0) side = ONSURF2OUT;
688 else side = ONSURF2IN;
689
690 return true;
691 }
692
693 /* ----------------------------------------------------------------------
694 check for axisymmetric move crossing line segment in (x,r) space
695 line segment from v1 to v2 in axisymmetry plane
696 3d move starting at x with v for tdelta
697 outface = exit face if particle hits no surface
698 used for special case test if vertical/horiz line segment is on exit face
699 solve quadratic eq to determine if curved trajectory intersects line seg
700 equivalent to straight trajectory intersecting cylindrical surf
701 can intersect infinite line 0,1,2 times
702 true intersection is at earliest positive time within line segment,
703 if all 3 conditions do not hold for 1st collision, can be 2nd collision
704 return 1 if yes, 0 if no
705 if yes, also return:
706 param = fraction of tdelta at collision pt
707 xc = x,y position of collision pt in axisymmetry plane
708 vc = vy,vz velocity at collision pt in axisymmetry plane
709 side = side of line that was hit = OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN
710 ------------------------------------------------------------------------- */
711
axi_line_intersect(double tdelta,double * x,double * v,int outface,double * lo,double * hi,double * v1,double * v2,double * norm,int selfflag,double * xc,double * vc,double & param,int & side)712 bool axi_line_intersect(double tdelta, double *x, double *v,
713 int outface, double *lo, double *hi,
714 double *v1, double *v2, double *norm, int selfflag,
715 double *xc, double *vc, double ¶m, int &side)
716 {
717 // compute nc = # of collisions with infinite line
718 // if 0, return false
719 // if 1, set t1
720 // if 2, set t1 and t2 with t1 < t2
721
722 int nc;
723 double t1,t2;
724
725 // vertical line segment
726 // no collision if starting on surface
727
728 if (v1[0] == v2[0]) {
729 nc = 1;
730 if (outface == 0 && v1[0] == lo[0]) t1 = tdelta;
731 else if (outface == 1 && v1[0] == hi[0]) t1 = tdelta;
732 else {
733 if (v[0] == 0.0) return false;
734 t1 = (v1[0] - x[0]) / v[0];
735 if (selfflag && t1 == 0.0) return false;
736 }
737
738 // horizontal line segment
739 // axi_horizontal_line() discards selfflag crossing
740
741 } else if (v1[1] == v2[1]) {
742 if (outface == 2 && v1[1] == lo[1]) {
743 nc = 1;
744 t1 = tdelta;
745 } else if (outface == 3 && v1[1] == hi[1]) {
746 nc = 1;
747 t1 = tdelta;
748 } else {
749 if (!axi_horizontal_line(tdelta,x,v,v1[1],nc,t1,t2)) return false;
750 }
751
752 // general line segment
753
754 } else {
755 double x21 = v2[0] - v1[0];
756 double y21 = v2[1] - v1[1];
757 double x21sq = x21*x21;
758 double y21sq = y21*y21;
759 double dconst = x21*v1[1] - y21*v1[0];
760
761 double a = x21sq*(v[1]*v[1] + v[2]*v[2]) - y21sq*v[0]*v[0];
762 if (a == 0.0) return false;
763 double b = x21sq*x[1]*v[1] - y21sq*x[0]*v[0] - y21*v[0]*dconst;
764 double c = x21sq*x[1]*x[1] - y21sq*x[0]*x[0] -
765 2.0*y21*x[0]*dconst - dconst*dconst;
766
767 double arg = b*b - a*c;
768 if (arg < 0.0) return false;
769 double sarg = sqrt(arg);
770
771 nc = 2;
772 double tone = (-b - sarg) / a;
773 double ttwo = (-b + sarg) / a;
774 t1 = MIN(tone,ttwo);
775 t2 = MAX(tone,ttwo);
776 }
777
778 // if selfflag, particle starts on surf line segment
779 // discard crossing at time = 0.0 or +/- epsilon (due to round-off)
780 // nc=2 test b/c horizontal line has already discarded this collision pt
781
782 if (selfflag && nc == 2) {
783 if (fabs(t1) < fabs(t2)) t1 = t2;
784 nc = 1;
785 }
786
787 // loop over 1 or 2 possible collision times
788
789 while (1) {
790
791 // test for collision time >= 0.0 and <= tdelta
792
793 if (t1 > tdelta) return false;
794 if (t1 < 0.0) {
795 if (nc == 1) return false;
796 t1 = t2;
797 nc--;
798 continue;
799 }
800
801 // NOTE: now doing this test above while loop, to avoid use of EPSSELF
802 // if selfflag, discard a collision near time = 0.0
803 // since is just a collision due to start on same surf just collided with
804 // reset selfflag = 0, so only do this for at most one of two collisions
805 // EPSSELF = 1.0e-6 seems to work, or 1.0e-4 and 1.0e-2
806 // 1.0e-8 does not unless push collision point to surf (below)
807
808 /*
809 if (selfflag && t1/tdelta < EPSSELF) {
810 if (nc == 1) return false;
811 selfflag = 0;
812 t1 = t2;
813 nc--;
814 continue;
815 }
816 */
817
818 // set xc[0,1] and vc[1,2] to values in axisymmetric plane
819 // if vertical or horizontal surf, insure xc is on it
820
821 xc[0] = x[0] + t1*v[0];
822 if (v1[0] == v2[0]) xc[0] = v1[0];
823 double ynew = x[1] + t1*v[1];
824 double znew = x[2] + t1*v[2];
825 xc[1] = sqrt(ynew*ynew + znew*znew);
826 if (v1[1] == v2[1]) xc[1] = v1[1];
827 xc[2] = 0.0;
828
829 double rn = ynew / xc[1];
830 double wn = znew / xc[1];
831 vc[0] = v[0];
832 vc[1] = v[1]*rn + v[2]*wn;
833 vc[2] = -v[1]*wn + v[2]*rn;
834
835 // test that xc is within line segment bounds
836 // y-test for vertical line, else x-test
837
838 bool within = true;
839 if (v1[0] == v2[0]) {
840 if (v1[1] < v2[1]) {
841 if (xc[1] < v1[1] || xc[1] > v2[1]) within = false;
842 } else {
843 if (xc[1] < v2[1] || xc[1] > v1[1]) within = false;
844 }
845 } else if (v1[0] < v2[0]) {
846 if (xc[0] < v1[0] || xc[0] > v2[0]) within = false;
847 } else {
848 if (xc[0] < v2[0] || xc[0] > v1[0]) within = false;
849 }
850
851 if (within) break;
852 if (nc == 1) return false;
853 t1 = t2;
854 nc--;
855 }
856
857 // could push collision point onto line segment
858 // so that selfflag test above could use smaller EPS
859 // could do this with projection operation instead
860
861 /*
862 double vec[3];
863 MathExtra::sub3(v2,v1,vec);
864 double lenbig = MathExtra::len3(vec);
865 MathExtra::sub3(xc,v1,vec);
866 double lensmall = MathExtra::len3(vec);
867 double ratio = lensmall/lenbig;
868 xc[0] = v1[0] + ratio*(v2[0]-v1[0]);
869 xc[1] = v1[1] + ratio*(v2[1]-v1[1]);
870 */
871
872 // there is a valid intersection with line segment
873 // set side to OUTSIDE or INSIDE
874 // no ONSURF case b/c surface is curved
875 // starting on surf is handled above by selfflag
876 // in axisymmetric plane, particle path is also curved
877 // regardless of where particle starts, it can hit front or back of surf
878 // use velocity vector at collision pt to determine side
879
880 double dot = MathExtra::dot3(norm,vc);
881 if (dot < 0.0) side = OUTSIDE;
882 else side = INSIDE;
883
884 param = t1/tdelta;
885 return true;
886 }
887
888 /* ----------------------------------------------------------------------
889 check for axisymmetric move crossing horizontal line in (x,r) space
890 not line segment but infinite horizontal line
891 called from Update for cell boundary
892 just uses first collision, must be within tdelta
893 can be called when particle is moving into cell from boundary
894 is not called when particle is moving out of cell from boundary (PEXIT)
895 called from axi_line_intersect() for special case of horizontal surf
896 uses one or two collisions, it will discard collisions outside segment
897 horizontal line is at yhoriz
898 move starting at x with v for tdelta
899 equations: solve for intersetion of line with circle in y,z plane
900 line: y = y0 + t*vy; z = t*vz since z0 = 0
901 circle: y^2 + z^2 = ylo^2
902 solve for t = time or intersection, via quadractic equation
903 special case in caller and here when y0 = ylo
904 otherwise 0 or 2 intersections of line with circle
905 first intersection has to happen within tdelta to be relevant
906 return 1 if crosses with 0.0 <= t <= tdelta
907 if crosses, also return nc, t1, t2 (can be two collisions)
908 ------------------------------------------------------------------------- */
909
axi_horizontal_line(double tdelta,double * x,double * v,double yhoriz,int & nc,double & t1,double & t2)910 bool axi_horizontal_line(double tdelta, double *x, double *v,
911 double yhoriz, int &nc, double &t1, double &t2)
912 {
913 double a = v[1]*v[1] + v[2]*v[2];
914 if (a == 0.0) return false;
915 double b = -v[1]*x[1];
916 double arg = yhoriz*yhoriz*a - v[2]*v[2]*x[1]*x[1];
917 if (arg < 0.0) return false;
918 double sarg = sqrt(arg);
919
920 nc = 2;
921 double tone = (b - sarg) / a;
922 double ttwo = (b + sarg) / a;
923 t1 = MIN(tone,ttwo);
924 t2 = MAX(tone,ttwo);
925
926 // if particle starts on line,
927 // discard crossing at time = 0.0 or +/- epsilon (due to round-off)
928 // due to cell crossing or selfflag in axi_line_intersect() caller
929
930 if (x[1] == yhoriz) {
931 if (fabs(t1) < fabs(t2)) t1 = t2;
932 nc = 1;
933 }
934
935 // require first collision time >= 0.0 and <= tdelta
936
937 if (t1 < 0.0 || t1 > tdelta) {
938 if (nc == 1) return false;
939 t1 = t2;
940 if (t1 < 0.0 || t1 > tdelta) return false;
941 nc = 1;
942 }
943
944 return true;
945 }
946
947 /* ----------------------------------------------------------------------
948 detect intersection between a directed line segment and a triangle
949 intersection is defined as any line segment pt (including end pts)
950 in common with any triangle pt (interior, edge, vertex)
951 one exception is if both line end pts are in plane of triangle,
952 then is NOT an intersection
953 start,stop = end points of directed line segment, can have zero length
954 v0,v1,v2 = 3 vertices of triangle
955 norm = unit vector normal to triangle plane
956 pointing OUTSIDE via right-hand rule
957 return TRUE if there is an intersection, else FALSE
958 if TRUE also return:
959 point = pt of intersection
960 param = intersection pt is this fraction along line (0-1 inclusive)
961 side = side of B that was hit = OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN
962 ------------------------------------------------------------------------- */
963
line_tri_intersect(double * start,double * stop,double * v0,double * v1,double * v2,double * norm,double * point,double & param,int & side)964 bool line_tri_intersect(double *start, double *stop,
965 double *v0, double *v1, double *v2, double *norm,
966 double *point, double ¶m, int &side)
967 {
968 double vec[3],start2stop[3],edge[3],pvec[3],xproduct[3];
969
970 // if start,stop are on same side of triangle, no intersection
971 // if start,stop are both in plane of triangle, no intersection
972
973 MathExtra::sub3(start,v0,vec);
974 double dotstart = MathExtra::dot3(norm,vec);
975 MathExtra::sub3(stop,v0,vec);
976 double dotstop = MathExtra::dot3(norm,vec);
977
978 if (dotstart < 0.0 && dotstop < 0.0) return false;
979 if (dotstart > 0.0 && dotstop > 0.0) return false;
980 if (dotstart == 0.0 && dotstop == 0.0) return false;
981
982 // param = parametric distance from start to stop
983 // at which tri plane is intersected
984 // force param to be 0.0 to 1.0 inclusive
985
986 MathExtra::sub3(v0,start,vec);
987 MathExtra::sub3(stop,start,start2stop);
988 param = MathExtra::dot3(norm,vec) / MathExtra::dot3(norm,start2stop);
989 param = MAX(param,0.0);
990 param = MIN(param,1.0);
991
992 // point = intersection pt with plane of triangle
993
994 point[0] = start[0] + param * start2stop[0];
995 point[1] = start[1] + param * start2stop[1];
996 point[2] = start[2] + param * start2stop[2];
997
998 // test if intersection pt is inside triangle
999 // edge = edge vector of triangle
1000 // pvec = vector from triangle vertex to intersection point
1001 // xproduct = cross product of edge with pvec
1002 // if dot product of xproduct with norm < 0.0 for any of 3 edges,
1003 // intersection point is outside tri
1004 // use EPSSQNEG instead of 0.0 for following case:
1005 // intersection pt is on tri edge where 2 tris come together
1006 // want it to detect collision with at least one of tris
1007 // point can be epsilon away from edge
1008 // this leads to xproduct being epsilon vec in opposite dirs for 2 tris
1009 // this can lead to dot3() being negative espilon^2 for both tris,
1010 // depending on direction of 2 tri norms
1011 // thus this can lead to no collision with either tri
1012 // typical observed dot values were -1.0e-18, so use EPSSQNEG = -1.0e-16
1013
1014 MathExtra::sub3(v1,v0,edge);
1015 MathExtra::sub3(point,v0,pvec);
1016 MathExtra::cross3(edge,pvec,xproduct);
1017 if (MathExtra::dot3(xproduct,norm) < EPSSQNEG) return false;
1018
1019 MathExtra::sub3(v2,v1,edge);
1020 MathExtra::sub3(point,v1,pvec);
1021 MathExtra::cross3(edge,pvec,xproduct);
1022 if (MathExtra::dot3(xproduct,norm) < EPSSQNEG) return false;
1023
1024 MathExtra::sub3(v0,v2,edge);
1025 MathExtra::sub3(point,v2,pvec);
1026 MathExtra::cross3(edge,pvec,xproduct);
1027 if (MathExtra::dot3(xproduct,norm) < EPSSQNEG) return false;
1028
1029 // there is a valid intersection with triangle
1030 // set side to ONSUFR, OUTSIDE, or INSIDE
1031 // if start point is inside or outside then side = same
1032 // if particle started on triangle, side = ONSURF OUT/IN based on dotstop
1033
1034 if (dotstart < 0.0) side = INSIDE;
1035 else if (dotstart > 0.0) side = OUTSIDE;
1036 else if (dotstop > 0.0) side = ONSURF2OUT;
1037 else side = ONSURF2IN;
1038
1039 return true;
1040 }
1041
1042 /* ----------------------------------------------------------------------
1043 determine which side of plane the point x,y,z is on
1044 plane is defined by vertex pt v and unit normal vec
1045 return -1,0,1 for below,on,above plane
1046 ------------------------------------------------------------------------- */
1047
whichside(double * v,double * norm,double x,double y,double z)1048 int whichside(double *v, double *norm, double x, double y, double z)
1049 {
1050 double vec[3];
1051 vec[0] = x - v[0];
1052 vec[1] = y - v[1];
1053 vec[2] = z - v[2];
1054
1055 double dotproduct = MathExtra::dot3(norm,vec);
1056 if (dotproduct < 0.0) return -1;
1057 else if (dotproduct > 0.0) return 1;
1058 else return 0;
1059 }
1060
1061 /* ----------------------------------------------------------------------
1062 determine if point x lies on surface of hex defined by lo and hi
1063 return 1 if it does, 0 if not
1064 ------------------------------------------------------------------------- */
1065
point_on_hex(double * x,double * lo,double * hi)1066 int point_on_hex(double *x, double *lo, double *hi)
1067 {
1068 if ((x[0] == lo[0] || x[0] == hi[0]) &&
1069 x[1] >= lo[1] && x[1] <= hi[1] && x[2] >= lo[2] && x[2] <= hi[2])
1070 return 1;
1071 if ((x[1] == lo[1] || x[1] == hi[1]) &&
1072 x[0] >= lo[0] && x[0] <= hi[0] && x[2] >= lo[2] && x[2] <= hi[2])
1073 return 1;
1074 if ((x[2] == lo[2] || x[2] == hi[2]) &&
1075 x[0] >= lo[0] && x[0] <= hi[0] && x[1] >= lo[1] && x[1] <= hi[1])
1076 return 1;
1077 return 0;
1078 }
1079
1080 /* ----------------------------------------------------------------------
1081 determine if point x lies inside or on surface of hex defined by lo and hi
1082 return 1 if it does, 0 if not
1083 ------------------------------------------------------------------------- */
1084
point_in_hex(double * x,double * lo,double * hi)1085 int point_in_hex(double *x, double *lo, double *hi)
1086 {
1087 if (x[0] >= lo[0] && x[0] <= hi[0] &&
1088 x[1] >= lo[1] && x[1] <= hi[1] &&
1089 x[2] >= lo[2] && x[2] <= hi[2]) return 1;
1090 return 0;
1091 }
1092
1093
1094 /* ----------------------------------------------------------------------
1095 determine if point x lies on edge or inside of tri defined by p1,p2,p3
1096 return 1 if it does, 0 if not
1097 ------------------------------------------------------------------------- */
1098
point_in_tri(double * x,double * p1,double * p2,double * p3,double * norm)1099 int point_in_tri(double *x, double *p1, double *p2, double *p3, double *norm)
1100 {
1101 // if not in plane of tri, then not inside tri
1102
1103 if (whichside(p1,norm,x[0],x[1],x[2])) return 0;
1104
1105 // enorm123 = 3 vecs normal to each edge of tri
1106 // are in plane of tri, pointing towards center of tri
1107 // enorms are NOT unit vectors
1108
1109 double enorm1[3],enorm2[3],enorm3[3];
1110
1111 double diff[3];
1112 MathExtra::sub3(p2,p1,diff);
1113 MathExtra::cross3(norm,diff,enorm1);
1114 MathExtra::sub3(p3,p2,diff);
1115 MathExtra::cross3(norm,diff,enorm2);
1116 MathExtra::sub3(p1,p3,diff);
1117 MathExtra::cross3(norm,diff,enorm3);
1118
1119 // if (pt - vertex) dotted into tri edge normal < 0, then outside tri
1120
1121 MathExtra::sub3(p1,x,diff);
1122 if (MathExtra::dot3(diff,enorm1) < 0.0) return 0;
1123 MathExtra::sub3(p2,x,diff);
1124 if (MathExtra::dot3(diff,enorm2) < 0.0) return 0;
1125 MathExtra::sub3(p3,x,diff);
1126 if (MathExtra::dot3(diff,enorm3) < 0.0) return 0;
1127 return 1;
1128 }
1129
1130 /* ----------------------------------------------------------------------
1131 compute distance bewteen a point X and line segment (P1,P2)
1132 distance = nearest distance to any point on line segment
1133 ------------------------------------------------------------------------- */
1134
distsq_point_line(double * x,double * p1,double * p2)1135 double distsq_point_line(double *x, double *p1, double *p2)
1136 {
1137 // A = vector from P1 to X
1138 // B = vector from P1 to P2
1139
1140 double a[3],b[3],c[3];
1141 MathExtra::sub3(x,p1,a);
1142 MathExtra::sub3(p2,p1,b);
1143
1144 // let P = projected point on infinite P1 to P2 line that is closest to X
1145 // alpha = fraction of distance from P1 to P2 that P is at
1146 // alpha can be < 0, or between 0 to 1, or > 1
1147
1148 double alpha = MathExtra::dot3(a,b)/MathExtra::lensq3(b);
1149
1150 // C = vector from point on P1P2 line to X
1151 // if alpha < 0.0, point on line is P1
1152 // if alpha > 1.0, point on line is P2
1153 // else point on line is P1 + alpha*(P2-P1)
1154
1155 if (alpha >= 1.0) MathExtra::sub3(x,p2,c);
1156 else if (alpha > 0.0) {
1157 a[0] = p1[0] + alpha*b[0];
1158 a[1] = p1[1] + alpha*b[1];
1159 a[2] = p1[2] + alpha*b[2];
1160 MathExtra::sub3(x,a,c);
1161 } else MathExtra::sub3(x,p1,c);
1162
1163 // return length of C
1164
1165 return MathExtra::lensq3(c);
1166 }
1167
1168 /* ----------------------------------------------------------------------
1169 compute distance bewteen a point X and triangle (P1,P2,P3) with NORM
1170 distance = nearest distance to any point within 2d triangle surface
1171 ------------------------------------------------------------------------- */
1172
distsq_point_tri(double * x,double * p1,double * p2,double * p3,double * norm)1173 double distsq_point_tri(double *x, double *p1, double *p2, double *p3,
1174 double *norm)
1175 {
1176 double a[3],point[3],edge[3],pvec[3],xproduct[3];
1177
1178 // A = vector from P1 to X
1179
1180 MathExtra::sub3(x,p1,a);
1181
1182 // point = projected point on infinite triangle plane
1183 // pdistsq = projected distance to plane
1184
1185 double alpha = MathExtra::dot3(a,norm);
1186 point[0] = x[0] - alpha*norm[0];
1187 point[1] = x[1] - alpha*norm[1];
1188 point[2] = x[2] - alpha*norm[2];
1189
1190 MathExtra::sub3(x,point,a);
1191 double pdistsq = MathExtra::lensq3(a);
1192
1193 // test if projected point is inside triangle
1194 // edge = edge vector of triangle
1195 // pvec = vector from triangle vertex to projected point
1196 // xproduct = cross product of edge with pvec
1197 // if dot product of xproduct with norm < 0.0 for any of 3 edges,
1198 // projected point is outside tri
1199 // if inside, return projected distance to plane
1200
1201 int inside = 1;
1202
1203 MathExtra::sub3(p2,p1,edge);
1204 MathExtra::sub3(point,p1,pvec);
1205 MathExtra::cross3(edge,pvec,xproduct);
1206 if (MathExtra::dot3(xproduct,norm) < 0.0) inside = 0;
1207
1208 MathExtra::sub3(p3,p2,edge);
1209 MathExtra::sub3(point,p2,pvec);
1210 MathExtra::cross3(edge,pvec,xproduct);
1211 if (MathExtra::dot3(xproduct,norm) < 0.0) inside = 0;
1212
1213 MathExtra::sub3(p1,p3,edge);
1214 MathExtra::sub3(point,p3,pvec);
1215 MathExtra::cross3(edge,pvec,xproduct);
1216 if (MathExtra::dot3(xproduct,norm) < 0.0) inside = 0;
1217
1218 if (inside) return pdistsq;
1219
1220 // projected point is outside triangle
1221 // compute minimum distance to any of 3 triangle edges
1222 // return sum of min distance and projected distance
1223
1224 double rsq = distsq_point_line(point,p1,p2);
1225 rsq = MIN(rsq,distsq_point_line(point,p2,p3));
1226 rsq = MIN(rsq,distsq_point_line(point,p3,p1));
1227 return rsq + pdistsq;
1228 }
1229
1230 /* ----------------------------------------------------------------------
1231 compute distance bewteen a line segmeht XY and 2d quad lo/hi
1232 ------------------------------------------------------------------------- */
1233
dist_line_quad(double * x,double * y,double * lo,double * hi)1234 double dist_line_quad(double *x, double *y, double *lo, double *hi)
1235 {
1236 double distsq;
1237 double pt[3],e1[3],e2[3];
1238
1239 pt[2] = e1[2] = e2[2] = 0.0;
1240
1241 // distance between 4 corner pts of quad and line segment
1242
1243 pt[0] = lo[0]; pt[1] = lo[1];
1244 distsq = distsq_point_line(pt,x,y);
1245 pt[0] = hi[0]; pt[1] = lo[1];
1246 distsq = MIN(distsq,distsq_point_line(pt,x,y));
1247 pt[0] = lo[0]; pt[1] = hi[1];
1248 distsq = MIN(distsq,distsq_point_line(pt,x,y));
1249 pt[0] = hi[0]; pt[1] = hi[1];
1250 distsq = MIN(distsq,distsq_point_line(pt,x,y));
1251
1252 // distance between line segment end pts and 4 quad edges
1253
1254 e1[0] = lo[0]; e1[1] = lo[1]; e2[0] = lo[0]; e2[1] = hi[1];
1255 distsq = MIN(distsq,distsq_point_line(x,e1,e2));
1256 distsq = MIN(distsq,distsq_point_line(y,e1,e2));
1257
1258 e1[0] = lo[0]; e1[1] = hi[1]; e2[0] = hi[0]; e2[1] = hi[1];
1259 distsq = MIN(distsq,distsq_point_line(x,e1,e2));
1260 distsq = MIN(distsq,distsq_point_line(y,e1,e2));
1261
1262 e1[0] = hi[0]; e1[1] = hi[1]; e2[0] = hi[0]; e2[1] = lo[1];
1263 distsq = MIN(distsq,distsq_point_line(x,e1,e2));
1264 distsq = MIN(distsq,distsq_point_line(y,e1,e2));
1265
1266 e1[0] = hi[0]; e1[1] = lo[1]; e2[0] = lo[0]; e2[1] = lo[1];
1267 distsq = MIN(distsq,distsq_point_line(x,e1,e2));
1268 distsq = MIN(distsq,distsq_point_line(y,e1,e2));
1269
1270 return sqrt(distsq);
1271 }
1272
1273 /* ----------------------------------------------------------------------
1274 compute distance bewteen a triangle XYZ with norm and 3d hex lo/hi
1275 ------------------------------------------------------------------------- */
1276
dist_tri_hex(double * x,double * y,double * z,double * norm,double * lo,double * hi)1277 double dist_tri_hex(double *x, double *y, double *z, double *norm,
1278 double *lo, double *hi)
1279 {
1280 double distsq;
1281 double pt[8][3],face[3];
1282
1283 // convert lo/hi to 8 corner pts
1284
1285 pt[0][0] = lo[0]; pt[0][1] = lo[1]; pt[0][2] = lo[2];
1286 pt[1][0] = hi[0]; pt[1][1] = lo[1]; pt[1][2] = lo[2];
1287 pt[2][0] = lo[0]; pt[2][1] = hi[1]; pt[2][2] = lo[2];
1288 pt[3][0] = hi[0]; pt[3][1] = hi[1]; pt[3][2] = lo[2];
1289 pt[4][0] = lo[0]; pt[4][1] = lo[1]; pt[4][2] = hi[2];
1290 pt[5][0] = hi[0]; pt[5][1] = lo[1]; pt[5][2] = hi[2];
1291 pt[6][0] = lo[0]; pt[6][1] = hi[1]; pt[6][2] = hi[2];
1292 pt[7][0] = hi[0]; pt[7][1] = hi[1]; pt[7][2] = hi[2];
1293
1294 // distance between 8 corner pts of hex and tri
1295
1296 distsq = distsq_point_tri(pt[0],x,y,z,norm);
1297 distsq = MIN(distsq,distsq_point_tri(pt[1],x,y,z,norm));
1298 distsq = MIN(distsq,distsq_point_tri(pt[2],x,y,z,norm));
1299 distsq = MIN(distsq,distsq_point_tri(pt[3],x,y,z,norm));
1300 distsq = MIN(distsq,distsq_point_tri(pt[4],x,y,z,norm));
1301 distsq = MIN(distsq,distsq_point_tri(pt[5],x,y,z,norm));
1302 distsq = MIN(distsq,distsq_point_tri(pt[6],x,y,z,norm));
1303 distsq = MIN(distsq,distsq_point_tri(pt[7],x,y,z,norm));
1304
1305 // distance between tri corner pts and 6 hex faces (2 tris per face)
1306
1307 face[0] = -1.0; face[1] = 0.0; face[2] = 0.0;
1308 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[4],pt[6],face));
1309 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[4],pt[6],face));
1310 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[4],pt[6],face));
1311 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[6],pt[2],face));
1312 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[6],pt[2],face));
1313 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[6],pt[2],face));
1314
1315 face[0] = 1.0; face[1] = 0.0; face[2] = 0.0;
1316 distsq = MIN(distsq,distsq_point_tri(x,pt[1],pt[3],pt[7],face));
1317 distsq = MIN(distsq,distsq_point_tri(y,pt[1],pt[3],pt[7],face));
1318 distsq = MIN(distsq,distsq_point_tri(z,pt[1],pt[3],pt[7],face));
1319 distsq = MIN(distsq,distsq_point_tri(x,pt[1],pt[7],pt[5],face));
1320 distsq = MIN(distsq,distsq_point_tri(y,pt[1],pt[7],pt[5],face));
1321 distsq = MIN(distsq,distsq_point_tri(z,pt[1],pt[7],pt[5],face));
1322
1323 face[0] = 0.0; face[1] = -1.0; face[2] = 0.0;
1324 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[1],pt[5],face));
1325 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[1],pt[5],face));
1326 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[1],pt[5],face));
1327 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[5],pt[4],face));
1328 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[5],pt[4],face));
1329 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[5],pt[4],face));
1330
1331 face[0] = 0.0; face[1] = 1.0; face[2] = 0.0;
1332 distsq = MIN(distsq,distsq_point_tri(x,pt[2],pt[6],pt[7],face));
1333 distsq = MIN(distsq,distsq_point_tri(y,pt[2],pt[6],pt[7],face));
1334 distsq = MIN(distsq,distsq_point_tri(z,pt[2],pt[6],pt[7],face));
1335 distsq = MIN(distsq,distsq_point_tri(x,pt[2],pt[7],pt[3],face));
1336 distsq = MIN(distsq,distsq_point_tri(y,pt[2],pt[7],pt[3],face));
1337 distsq = MIN(distsq,distsq_point_tri(z,pt[2],pt[7],pt[3],face));
1338
1339 face[0] = 0.0; face[1] = 0.0; face[2] = -1.0;
1340 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[2],pt[3],face));
1341 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[2],pt[3],face));
1342 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[2],pt[3],face));
1343 distsq = MIN(distsq,distsq_point_tri(x,pt[0],pt[3],pt[1],face));
1344 distsq = MIN(distsq,distsq_point_tri(y,pt[0],pt[3],pt[1],face));
1345 distsq = MIN(distsq,distsq_point_tri(z,pt[0],pt[3],pt[1],face));
1346
1347 face[0] = 0.0; face[1] = 0.0; face[2] = 1.0;
1348 distsq = MIN(distsq,distsq_point_tri(x,pt[4],pt[5],pt[7],face));
1349 distsq = MIN(distsq,distsq_point_tri(y,pt[4],pt[5],pt[7],face));
1350 distsq = MIN(distsq,distsq_point_tri(z,pt[4],pt[5],pt[7],face));
1351 distsq = MIN(distsq,distsq_point_tri(x,pt[4],pt[7],pt[6],face));
1352 distsq = MIN(distsq,distsq_point_tri(y,pt[4],pt[7],pt[6],face));
1353 distsq = MIN(distsq,distsq_point_tri(z,pt[4],pt[7],pt[6],face));
1354
1355 return sqrt(distsq);
1356 }
1357
1358 /* ----------------------------------------------------------------------
1359 return minimum fractional distance that X is from line pts V0 or V1
1360 lensq = length of line
1361 frac = (length of X-V0 or X-V1) / lensq
1362 return fracsq = square of frac
1363 ------------------------------------------------------------------------- */
1364
line_fraction(double * x,double * v0,double * v1)1365 double line_fraction(double *x, double *v0, double *v1)
1366 {
1367 double segment[3];
1368
1369 MathExtra::sub3(v0,v1,segment);
1370 double lensq = MathExtra::lensq3(segment);
1371
1372 MathExtra::sub3(x,v0,segment);
1373 double fracsq = MathExtra::lensq3(segment)/lensq;
1374 MathExtra::sub3(x,v1,segment);
1375 fracsq = MIN(fracsq,MathExtra::lensq3(segment)/lensq);
1376
1377 return fracsq;
1378 }
1379
1380 /* ----------------------------------------------------------------------
1381 return minimum fractional distance that X is from triangle pts V0, V1, V2
1382 lensq = min length of any of 3 triangle edges
1383 frac = (length of X-V0 or X-V1 or X-V2) / lensq
1384 return fracsq = square of frac
1385 ------------------------------------------------------------------------- */
1386
tri_fraction(double * x,double * v0,double * v1,double * v2)1387 double tri_fraction(double *x, double *v0, double *v1, double *v2)
1388 {
1389 double segment[3];
1390
1391 MathExtra::sub3(v0,v1,segment);
1392 double lensq = MathExtra::lensq3(segment);
1393 MathExtra::sub3(v1,v2,segment);
1394 lensq = MIN(lensq,MathExtra::lensq3(segment));
1395 MathExtra::sub3(v0,v2,segment);
1396 lensq = MIN(lensq,MathExtra::lensq3(segment));
1397
1398 MathExtra::sub3(x,v0,segment);
1399 double fracsq = MathExtra::lensq3(segment)/lensq;
1400 MathExtra::sub3(x,v1,segment);
1401 fracsq = MIN(fracsq,MathExtra::lensq3(segment)/lensq);
1402 MathExtra::sub3(x,v2,segment);
1403 fracsq = MIN(fracsq,MathExtra::lensq3(segment)/lensq);
1404
1405 return fracsq;
1406 }
1407
1408 /* ---------------------------------------------------------------------- */
1409
1410 }
1411