Block-Structured AMR Software Framework
 
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
AMReX_EB_triGeomOps_K.H
Go to the documentation of this file.
1#ifndef AMREX_EB_TRIGEOMOPS_K_H_
2#define AMREX_EB_TRIGEOMOPS_K_H_
3#include <AMReX_Config.H>
4#include <AMReX_Math.H>
5#include <AMReX_REAL.H>
6
7namespace amrex::tri_geom_ops
8{
9 //================================================================================
10 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real Distance2(const Real P1[3],const Real P2[3])
11 {
12 return( (P1[0]-P2[0])*(P1[0]-P2[0]) +
13 (P1[1]-P2[1])*(P1[1]-P2[1]) +
14 (P1[2]-P2[2])*(P1[2]-P2[2]) );
15 }
16 //================================================================================
17 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real DotProd(const Real v1[3],const Real v2[3])
18 {
19 return(v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
20 }
21 //================================================================================
22 // If this is zero, the two lines either intersect or are parallel.
23 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real side_op(const Real L1[6],const Real L2[6])
24 {
25 return( L1[0]*L2[4]
26 + L1[1]*L2[5]
27 + L1[2]*L2[3]
28 + L1[3]*L2[2]
29 + L1[4]*L2[0]
30 + L1[5]*L2[1] );
31 }
32 //================================================================================
33 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getvec(const Real P1[3],const Real P2[3],Real v[3])
34 {
35 v[0]=P2[0]-P1[0];
36 v[1]=P2[1]-P1[1];
37 v[2]=P2[2]-P1[2];
38 }
39 //================================================================================
40 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getunitvec(const Real v[3],Real vu[3])
41 {
42 Real vmag;
43 vmag=std::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
44 vu[0]=v[0]/vmag;
45 vu[1]=v[1]/vmag;
46 vu[2]=v[2]/vmag;
47 }
48 //================================================================================
49 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void CrossProd(const Real v1[3],const Real v2[3],Real v[3])
50 {
51 v[0]=v1[1]*v2[2]-v1[2]*v2[1];
52 v[1]=v1[2]*v2[0]-v1[0]*v2[2];
53 v[2]=v1[0]*v2[1]-v1[1]*v2[0];
54 }
55 //================================================================================
56 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void get_plucker_coords(const Real v1[3],const Real v2[3],Real L[6])
57 {
58 L[0] = v1[0]*v2[1] - v1[1]*v2[0];
59 L[1] = v1[0]*v2[2] - v1[2]*v2[0];
60 L[2] = v1[0] - v2[0];
61 L[3] = v1[1]*v2[2] - v1[2]*v2[1];
62 L[4] = v1[2] - v2[2];
63 L[5] = v2[1] - v1[1];
64 }
65 //================================================================================
66 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void side_op3(const Real v1[3], const Real v2[3],
67 const Real t1[3], const Real t2[3], const Real t3[3],
68 Real &S1, Real &S2, Real &S3)
69 {
70
71 Real L[6],e1[6],e2[6],e3[6];
72
73 get_plucker_coords(v1,v2,L);
74 get_plucker_coords(t1,t2,e1);
75 get_plucker_coords(t2,t3,e2);
76 get_plucker_coords(t3,t1,e3);
77
78 S1=side_op(L,e1);
79 S2=side_op(L,e2);
80 S3=side_op(L,e3);
81 }
82 //================================================================================
83 //get normal of triangle pointing at a test-point
84 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void tri_n(const Real P1[3], const Real P2[3], const Real P3[3],
85 const Real testp[3], Real n[3])
86 {
87 Real v1[3],v2[3],magn;
88 Real centr[3],c_tp_vec[3];
89
90 getvec(P1,P2,v1);
91 getvec(P1,P3,v2);
92 CrossProd(v1,v2,n);
93
94
95 centr[0]=Real(1./3.)*(P1[0]+P2[0]+P3[0]);
96 centr[1]=Real(1./3.)*(P1[1]+P2[1]+P3[1]);
97 centr[2]=Real(1./3.)*(P1[2]+P2[2]+P3[2]);
98
99 getvec(centr,testp,c_tp_vec);
100 magn=std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
101
102 if(DotProd(c_tp_vec,n) < Real(0.0))
103 {
104 magn=-magn;
105 }
106
107 n[0]=n[0]/magn;
108 n[1]=n[1]/magn;
109 n[2]=n[2]/magn;
110 }
111 //================================================================================
112 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real triangle_area(const Real P1[3], const Real P2[3], const Real P3[3])
113 {
114 Real v1[3],v2[3],area[3];
115
116 getvec(P1,P2,v1);
117 getvec(P1,P3,v2);
118 CrossProd(v1,v2,area);
119 return(Real(0.5) * std::sqrt(DotProd(area,area)) );
120 }
121 //================================================================================
122 //this is only useful when v1-v2 segment intersects the triangle
123 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool find_intersection_point(const Real v1[3],const Real v2[3],
124 const Real t1[3], const Real t2[3], const Real t3[3], Real ip[3],int bisect_iters=20,Real tol=1e-6)
125 {
126 Real plane_eq_mid,plane_eq1,plane_eq2;
127
128 Real ab[3],ac[3],n[3],magn;
129 Real midp[3],p1[3],p2[3];
130
131 getvec(t1,t2,ab);
132 getvec(t1,t3,ac);
133
134 CrossProd(ab,ac,n);
135 magn=std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
136
137 n[0]=n[0]/magn;
138 n[1]=n[1]/magn;
139 n[2]=n[2]/magn;
140
141 p1[0]=v1[0];
142 p1[1]=v1[1];
143 p1[2]=v1[2];
144
145 p2[0]=v2[0];
146 p2[1]=v2[1];
147 p2[2]=v2[2];
148
149
150 bool all_ok=true;
151
152 for(int i=0;i<bisect_iters;i++)
153 {
154 midp[0]=Real(0.5)*(p1[0]+p2[0]);
155 midp[1]=Real(0.5)*(p1[1]+p2[1]);
156 midp[2]=Real(0.5)*(p1[2]+p2[2]);
157
158 plane_eq_mid= (midp[0]-t1[0])*n[0] + (midp[1]-t1[1])*n[1] + (midp[2]-t1[2])*n[2];
159 plane_eq1 = (p1[0] -t1[0])*n[0] + (p1[1] -t1[1])*n[1] + (p1[2] -t1[2])*n[2];
160 plane_eq2 = (p2[0] -t1[0])*n[0] + (p2[1] -t1[1])*n[1] + (p2[2] -t1[2])*n[2];
161
162 //Print()<<"midp:"<<midp[0]<<"\t"<<midp[1]<<"\t"<<midp[2]<<"\t"<<plane_eq_mid<<"\n";
163
164 if(std::abs(plane_eq_mid) < tol)
165 {
166 break;
167 }
168
169 if(plane_eq_mid*plane_eq1 < Real(0.0))
170 {
171 p2[0]=midp[0];
172 p2[1]=midp[1];
173 p2[2]=midp[2];
174 }
175 else if(plane_eq_mid*plane_eq2 < Real(0.0))
176 {
177 p1[0]=midp[0];
178 p1[1]=midp[1];
179 p1[2]=midp[2];
180 }
181 else //plane_eq_mid is 0
182 //or error: p1,midp and p2 are on the same side
183 //which is not what this function is meant for
184 {
185 if(plane_eq_mid*plane_eq1 > Real(0.0) && plane_eq_mid*plane_eq2 > Real(0.0))
186 {
187 all_ok=false;
188 }
189 break;
190 }
191 }
192
193 ip[0]=midp[0];
194 ip[1]=midp[1];
195 ip[2]=midp[2];
196
197 return(all_ok);
198 }
199 //================================================================================
200 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int lineseg_tri_intersect(const Real v1[3], const Real v2[3],
201 const Real t1[3], const Real t2[3], const Real t3[3])
202 {
203 //see plucker coordinates based method
204 //https://members.loria.fr/SLazard/ARC-Visi3D/Pant-project/files/Line_Triangle.html
205
206 Real S1,S2,S3;
207 Real tri_area,area1,area2;
208 Real L2[6],L3[6],L4[6],ls_s1,ls_s2;
209
210 side_op3(v1,v2,t1,t2,t3,S1,S2,S3);
211
212 //we are assuming there are no intersections initially
213 int no_intersections=1;
214
215 Real eps = std::numeric_limits<Real>::epsilon();
216
217 //coplanar (S1,S2,S3 = 0)
218 if(std::abs(S1) < eps && std::abs(S2) < eps && std::abs(S3) < eps)
219 {
220 //Print()<<"line segment and triangle are in the same plane\t"<<S1<<"\t"<<S2<<"\t"<<S3<<"\n";
221 tri_area=triangle_area(t1,t2,t3);
222
223 /*if(tri_area == 0)
224 {
225 amrex::Abort("problem with triangle\n");
226 }*/
227 area1=(triangle_area(t1,t2,v1)+triangle_area(t2,t3,v1)+triangle_area(t3,t1,v1));
228 area2=(triangle_area(t1,t2,v2)+triangle_area(t2,t3,v2)+triangle_area(t3,t1,v2));
229
230 if( std::abs(area1-tri_area)>eps || std::abs(area2-tri_area)>eps )
231 {
232 no_intersections = 0;
233 }
234 }
235 //proper and edge intersection
236 else if( (S1 < Real(0.0) && S2 < Real(0.0) && S3 < Real(0.0)) ||
237 (S1 > Real(0.0) && S2 > Real(0.0) && S3 > Real(0.0)) ||
238 (std::abs(S1) < eps && S2*S3 > Real(0.0)) || //S1=0
239 (std::abs(S2) < eps && S3*S1 > Real(0.0)) || //S2=0
240 (std::abs(S3) < eps && S1*S2 > Real(0.0)) ) //S3=0
241 {
242
243 get_plucker_coords(v1,t1,L2);
244 get_plucker_coords(t1,v2,L3);
245 get_plucker_coords(t2,t3,L4);
246
247 /*if(std::abs(S1*S2*S3) < eps)
248 {
249 Print()<<"edge intersection S1,S2,S3:"
250 <<S1<<"\t"<<S2<<"\t"<<S3<<"\n";
251 }*/
252
253 ls_s1 = side_op(L4,L3);
254 ls_s2 = side_op(L4,L2);
255
256 if(ls_s1*ls_s2 > Real(0.0))
257 {
258 no_intersections = 0;
259 }
260 }
261 //vertex intersection
262 else if( (std::abs(S1) < eps && std::abs(S2) < eps) || //S1,S2=0
263 (std::abs(S2) < eps && std::abs(S3) < eps) ) //S2,S3=0
264 {
265
266 //Print()<<"vertex intersection type 1\n";
267 //don't chose vertex 2 or 3
268 get_plucker_coords(v2,t1,L2);
269 get_plucker_coords(t1,v1,L3);
270 get_plucker_coords(t2,t3,L4);
271
272 ls_s1=side_op(L4,L3);
273 ls_s2=side_op(L4,L2);
274
275 if(ls_s1*ls_s2 > 0)
276 {
277 no_intersections = 0;
278 }
279 }
280 else if(std::abs(S3) < eps && std::abs(S1) < eps) //S3,S1=0
281 {
282
283 //Print()<<"vertex intersection type 2\n";
284 //don't chose vertex 1
285 get_plucker_coords(v2,t2,L2);
286 get_plucker_coords(t2,v1,L3);
287 get_plucker_coords(t3,t1,L4);
288
289 ls_s1=side_op(L4,L3);
290 ls_s2=side_op(L4,L2);
291
292 if(ls_s1*ls_s2 > 0)
293 {
294 no_intersections=0;
295 }
296 }
297
298 return(no_intersections);
299
300 }
301 //================================================================================
302}
303#endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
AMReX_Math.H
AMReX_REAL.H
amrex::tri_geom_ops
Definition AMReX_EB_triGeomOps_K.H:8
amrex::tri_geom_ops::side_op
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real side_op(const Real L1[6], const Real L2[6])
Definition AMReX_EB_triGeomOps_K.H:23
amrex::tri_geom_ops::DotProd
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real DotProd(const Real v1[3], const Real v2[3])
Definition AMReX_EB_triGeomOps_K.H:17
amrex::tri_geom_ops::lineseg_tri_intersect
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int lineseg_tri_intersect(const Real v1[3], const Real v2[3], const Real t1[3], const Real t2[3], const Real t3[3])
Definition AMReX_EB_triGeomOps_K.H:200
amrex::tri_geom_ops::CrossProd
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void CrossProd(const Real v1[3], const Real v2[3], Real v[3])
Definition AMReX_EB_triGeomOps_K.H:49
amrex::tri_geom_ops::get_plucker_coords
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void get_plucker_coords(const Real v1[3], const Real v2[3], Real L[6])
Definition AMReX_EB_triGeomOps_K.H:56
amrex::tri_geom_ops::getvec
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getvec(const Real P1[3], const Real P2[3], Real v[3])
Definition AMReX_EB_triGeomOps_K.H:33
amrex::tri_geom_ops::side_op3
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void side_op3(const Real v1[3], const Real v2[3], const Real t1[3], const Real t2[3], const Real t3[3], Real &S1, Real &S2, Real &S3)
Definition AMReX_EB_triGeomOps_K.H:66
amrex::tri_geom_ops::Distance2
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real Distance2(const Real P1[3], const Real P2[3])
Definition AMReX_EB_triGeomOps_K.H:10
amrex::tri_geom_ops::getunitvec
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void getunitvec(const Real v[3], Real vu[3])
Definition AMReX_EB_triGeomOps_K.H:40
amrex::tri_geom_ops::tri_n
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void tri_n(const Real P1[3], const Real P2[3], const Real P3[3], const Real testp[3], Real n[3])
Definition AMReX_EB_triGeomOps_K.H:84
amrex::tri_geom_ops::triangle_area
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real triangle_area(const Real P1[3], const Real P2[3], const Real P3[3])
Definition AMReX_EB_triGeomOps_K.H:112
amrex::tri_geom_ops::find_intersection_point
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool find_intersection_point(const Real v1[3], const Real v2[3], const Real t1[3], const Real t2[3], const Real t3[3], Real ip[3], int bisect_iters=20, Real tol=1e-6)
Definition AMReX_EB_triGeomOps_K.H:123