Skip to content

Commit 65bf80d

Browse files
WeiqunZhangWeiqunZhang
authored
Add STLtools::fillSignedDistance (#4501)
1 parent 335f529 commit 65bf80d

File tree

4 files changed

+292
-2
lines changed

4 files changed

+292
-2
lines changed

‎Src/Base/AMReX_Dim3.H‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,30 @@ namespace amrex {
1212
struct Dim3 { int x; int y; int z; };
1313
struct XDim3 { Real x; Real y; Real z; };
1414

15+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
16+
Real dot_product (XDim3 const& a, XDim3 const& b)
17+
{
18+
return a.x*b.x + a.y*b.y + a.z*b.z;
19+
}
20+
21+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
22+
XDim3 cross_product (XDim3 const& a, XDim3 const& b)
23+
{
24+
return {a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x};
25+
}
26+
27+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
28+
XDim3 operator+ (XDim3 const& a, XDim3 const& b)
29+
{
30+
return {a.x + b.x, a.y + b.y, a.z + b.z};
31+
}
32+
33+
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
34+
XDim3 operator- (XDim3 const& a, XDim3 const& b)
35+
{
36+
return {a.x - b.x, a.y - b.y, a.z - b.z};
37+
}
38+
1539
template <typename T,
1640
std::enable_if_t <std::is_same_v<T, Dim3> ||
1741
std::is_same_v<T, XDim3>>* = nullptr>

‎Src/EB/AMReX_EBToPVD.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
#include <algorithm>
1010
#include <limits>
1111

12+
namespace amrex {
13+
1214
namespace {
1315
amrex::Real dot_product(const std::array<amrex::Real,3>& a, const std::array<amrex::Real,3>& b)
1416
{
@@ -22,8 +24,6 @@ bool intersects(amrex::Real val)
2224

2325
}
2426

25-
namespace amrex {
26-
2727
void EBToPVD::EBToPolygon(const Real* problo, const Real* dx,
2828
const Box & bx, Array4<EBCellFlag const> const& flag,
2929
Array4<Real const> const& bcent,

‎Src/EB/AMReX_EB_STL_utils.H‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,8 @@ public:
8484
Array<Array4<EB2::Type_t const>,AMREX_SPACEDIM> const& type_arr,
8585
Array4<Real const> const& lst, Geometry const& geom) ;
8686

87+
void fillSignedDistance (MultiFab& mf, IntVect const& nghost, Geometry const& geom) const;
88+
8789
void prepare (Gpu::PinnedVector<Triangle> a_tri_pts); // public for cuda
8890

8991
private:

‎Src/EB/AMReX_EB_STL_utils.cpp‎

Lines changed: 264 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,217 @@ namespace {
189189
}
190190
}
191191
}
192+
193+
#if (AMREX_SPACEDIM == 3)
194+
195+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
196+
Real pt_segment_min_d2 (XDim3 const& pt, XDim3 const& a, XDim3 const& b)
197+
{
198+
auto ab = b - a;
199+
auto ap = pt - a;
200+
auto t = dot_product(ab,ap) / dot_product(ab,ab);
201+
t = amrex::Clamp(t, Real(0), Real(1));
202+
return Math::powi<2>(t*ab.x-ap.x)
203+
+ Math::powi<2>(t*ab.y-ap.y)
204+
+ Math::powi<2>(t*ab.z-ap.z);
205+
}
206+
207+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
208+
Real pt_tri_min_d2 (XDim3 const& pt, STLtools::Triangle const& tri)
209+
{
210+
// We want to know whether the projection of the point onto the
211+
// triangle plane is inside the triangle. If so, the closest point
212+
// on the triangle is that projected point. If not, the closest
213+
// point must be on the edges or vertices of the triangle.
214+
215+
// Any point Q on the plane of a triangle (ABC) can be written as Q
216+
// = c1*A + c2*B + c3*C, where c1 + c2 + c3 = 1. The coefficients,
217+
// c1, c2 and c3, are also known as barycentric coordinates of point
218+
// Q relative to the triangle. If all coefficients are between 0 and
219+
// 1, the point lies inisde the triangle. The way to see this is
220+
// that point Q is the interpolation of the three vertices in that
221+
// case. Any point Q on the plance can also be written as (Q-A) =
222+
// s*(B-A) + t*(C-A). The relationship between the two sets of
223+
// coefficients are c1 = 1-s-t, c2 = s, and c3 = t.
224+
225+
// Let AQ = s*AB + t*AC. Let AP = AQ + d*n/|n|, where n is a vector
226+
// perpendicular to the plane, and d is the distance to the
227+
// plane. We can then solve a system of equations.
228+
auto vb = tri.v2 - tri.v1; // AB, i.e., B-A
229+
auto vc = tri.v3 - tri.v1; // AC
230+
auto n = cross_product(vb, vc);
231+
auto rhs = pt - tri.v1; // AP
232+
auto det = dot_product(n,n);
233+
auto detinv = Real(1) / det;
234+
auto s = dot_product(rhs, cross_product(vc, n)) * detinv;
235+
if (s >= 0 && s <= Real(1)) {
236+
auto t = dot_product(vb, cross_product(rhs, n)) * detinv;
237+
auto c1 = Real(1) - s - t;
238+
if (t >= 0 && c1 >= 0) { // q is inside
239+
auto dn = dot_product(vb, cross_product(vc, rhs));
240+
return dn*dn*detinv;
241+
}
242+
}
243+
244+
auto d2_1 = pt_segment_min_d2(pt, tri.v1, tri.v2);
245+
auto d2_2 = pt_segment_min_d2(pt, tri.v1, tri.v3);
246+
auto d2_3 = pt_segment_min_d2(pt, tri.v2, tri.v3);
247+
return std::min({d2_1,d2_2,d2_3});
248+
}
249+
250+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
251+
Real pt_box_min_d2 (XDim3 const& pt, RealBox const& bbox)
252+
{
253+
bool inside_x = (pt.x >= bbox.lo(0)) && (pt.x <= bbox.hi(0));
254+
bool inside_y = (pt.y >= bbox.lo(1)) && (pt.y <= bbox.hi(1));
255+
bool inside_z = (pt.z >= bbox.lo(2)) && (pt.z <= bbox.hi(2));
256+
if (inside_x && inside_y && inside_z) {
257+
return Real(0);
258+
} else if (inside_x && inside_y) {
259+
if (pt.z < bbox.lo(2)) {
260+
return Math::powi<2>(pt.z-bbox.lo(2));
261+
} else {
262+
return Math::powi<2>(pt.z-bbox.hi(2));
263+
}
264+
} else if (inside_x && inside_z) {
265+
if (pt.y < bbox.lo(1)) {
266+
return Math::powi<2>(pt.y-bbox.lo(1));
267+
} else {
268+
return Math::powi<2>(pt.y-bbox.hi(1));
269+
}
270+
} else if (inside_y && inside_z) {
271+
if (pt.x < bbox.lo(0)) {
272+
return Math::powi<2>(pt.x-bbox.lo(0));
273+
} else {
274+
return Math::powi<2>(pt.x-bbox.hi(0));
275+
}
276+
} else if (inside_x) {
277+
if ((pt.y < bbox.lo(1)) && (pt.z < bbox.lo(2))) {
278+
return Math::powi<2>(pt.y-bbox.lo(1))
279+
+ Math::powi<2>(pt.z-bbox.lo(2));
280+
} else if ((pt.y < bbox.lo(1)) && (pt.z > bbox.hi(2))) {
281+
return Math::powi<2>(pt.y-bbox.lo(1))
282+
+ Math::powi<2>(pt.z-bbox.hi(2));
283+
} else if ((pt.y > bbox.hi(1)) && (pt.z < bbox.lo(2))) {
284+
return Math::powi<2>(pt.y-bbox.hi(1))
285+
+ Math::powi<2>(pt.z-bbox.lo(2));
286+
} else {
287+
return Math::powi<2>(pt.y-bbox.hi(1))
288+
+ Math::powi<2>(pt.z-bbox.hi(2));
289+
}
290+
} else if (inside_y) {
291+
if ((pt.x < bbox.lo(0)) && (pt.z < bbox.lo(2))) {
292+
return Math::powi<2>(pt.x-bbox.lo(0))
293+
+ Math::powi<2>(pt.z-bbox.lo(2));
294+
} else if ((pt.x < bbox.lo(0)) && (pt.z > bbox.hi(2))) {
295+
return Math::powi<2>(pt.x-bbox.lo(0))
296+
+ Math::powi<2>(pt.z-bbox.hi(2));
297+
} else if ((pt.x > bbox.hi(0)) && (pt.z < bbox.lo(2))) {
298+
return Math::powi<2>(pt.x-bbox.hi(0))
299+
+ Math::powi<2>(pt.z-bbox.lo(2));
300+
} else {
301+
return Math::powi<2>(pt.x-bbox.hi(0))
302+
+ Math::powi<2>(pt.z-bbox.hi(2));
303+
}
304+
} else if (inside_z) {
305+
if ((pt.x < bbox.lo(0)) && (pt.y < bbox.lo(1))) {
306+
return Math::powi<2>(pt.x-bbox.lo(0))
307+
+ Math::powi<2>(pt.y-bbox.lo(1));
308+
} else if ((pt.x < bbox.lo(0)) && (pt.y > bbox.hi(1))) {
309+
return Math::powi<2>(pt.x-bbox.lo(0))
310+
+ Math::powi<2>(pt.y-bbox.hi(1));
311+
312+
} else if ((pt.x > bbox.hi(0)) && (pt.y < bbox.lo(1))) {
313+
return Math::powi<2>(pt.x-bbox.hi(0))
314+
+ Math::powi<2>(pt.y-bbox.lo(1));
315+
316+
} else {
317+
return Math::powi<2>(pt.x-bbox.hi(0))
318+
+ Math::powi<2>(pt.y-bbox.hi(1));
319+
320+
}
321+
} else {
322+
if ((pt.x < bbox.lo(0)) && (pt.y < bbox.lo(1)) && (pt.z < bbox.lo(2))) {
323+
return Math::powi<2>(pt.x-bbox.lo(0))
324+
+ Math::powi<2>(pt.y-bbox.lo(1))
325+
+ Math::powi<2>(pt.z-bbox.lo(2));
326+
} else if ((pt.x > bbox.hi(0)) && (pt.y < bbox.lo(1)) && (pt.z < bbox.lo(2))) {
327+
return Math::powi<2>(pt.x-bbox.hi(0))
328+
+ Math::powi<2>(pt.y-bbox.lo(1))
329+
+ Math::powi<2>(pt.z-bbox.lo(2));
330+
} else if ((pt.x < bbox.lo(0)) && (pt.y > bbox.hi(1)) && (pt.z < bbox.lo(2))) {
331+
return Math::powi<2>(pt.x-bbox.lo(0))
332+
+ Math::powi<2>(pt.y-bbox.hi(1))
333+
+ Math::powi<2>(pt.z-bbox.lo(2));
334+
} else if ((pt.x > bbox.hi(0)) && (pt.y > bbox.hi(1)) && (pt.z < bbox.lo(2))) {
335+
return Math::powi<2>(pt.x-bbox.hi(0))
336+
+ Math::powi<2>(pt.y-bbox.hi(1))
337+
+ Math::powi<2>(pt.z-bbox.lo(2));
338+
} else if ((pt.x < bbox.lo(0)) && (pt.y < bbox.lo(1)) && (pt.z > bbox.hi(2))) {
339+
return Math::powi<2>(pt.x-bbox.lo(0))
340+
+ Math::powi<2>(pt.y-bbox.lo(1))
341+
+ Math::powi<2>(pt.z-bbox.hi(2));
342+
} else if ((pt.x > bbox.hi(0)) && (pt.y < bbox.lo(1)) && (pt.z > bbox.hi(2))) {
343+
return Math::powi<2>(pt.x-bbox.hi(0))
344+
+ Math::powi<2>(pt.y-bbox.lo(1))
345+
+ Math::powi<2>(pt.z-bbox.hi(2));
346+
} else if ((pt.x < bbox.lo(0)) && (pt.y > bbox.hi(1)) && (pt.z > bbox.hi(2))) {
347+
return Math::powi<2>(pt.x-bbox.lo(0))
348+
+ Math::powi<2>(pt.y-bbox.hi(1))
349+
+ Math::powi<2>(pt.z-bbox.hi(2));
350+
} else {
351+
return Math::powi<2>(pt.x-bbox.hi(0))
352+
+ Math::powi<2>(pt.y-bbox.hi(1))
353+
+ Math::powi<2>(pt.z-bbox.hi(2));
354+
}
355+
}
356+
}
357+
358+
template <int M, int N>
359+
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
360+
Real bvh_d2 (XDim3 const& pt, STLtools::BVHNodeT<M,N> const* root)
361+
{
362+
Stack<int, STLtools::m_bvh_max_stack_size> nodes_to_do;
363+
Stack<std::int8_t, STLtools::m_bvh_max_stack_size> nchildren_done;
364+
nodes_to_do.push(0);
365+
nchildren_done.push(0);
366+
367+
Real d = std::numeric_limits<Real>::max();
368+
369+
while (!nodes_to_do.empty()) {
370+
auto const& node = root[nodes_to_do.top()];
371+
if (node.nchildren == 0) { // leaf node
372+
for (std::int8_t it = 0; it < node.ntriangles; ++it) {
373+
Real dmin = pt_tri_min_d2(pt, node.triangles[it]);
374+
d = std::min(d,dmin);
375+
if (d == 0) { return 0; }
376+
}
377+
nodes_to_do.pop();
378+
nchildren_done.pop();
379+
} else {
380+
auto& ndone = nchildren_done.top();
381+
if (ndone < node.nchildren) {
382+
for (auto ichild = ndone; ichild < node.nchildren; ++ichild) {
383+
++ndone;
384+
int inode = node.children[ichild];
385+
auto dmin = pt_box_min_d2(pt, root[inode].boundingbox);
386+
if (d > dmin) {
387+
nodes_to_do.push(inode);
388+
nchildren_done.push(0);
389+
break;
390+
}
391+
}
392+
} else {
393+
nodes_to_do.pop();
394+
nchildren_done.pop();
395+
}
396+
}
397+
}
398+
399+
return d;
400+
}
401+
402+
#endif
192403
}
193404

194405
void
@@ -1113,4 +1324,57 @@ STLtools::updateIntercept (Array<Array4<Real>,AMREX_SPACEDIM> const& inter_arr,
11131324
}
11141325
}
11151326

1327+
void
1328+
STLtools::fillSignedDistance (MultiFab& mf, IntVect const& nghost, Geometry const& geom) const
1329+
{
1330+
BL_PROFILE("STLtools::fillSignedDistance");
1331+
1332+
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(AMREX_SPACEDIM == 3,
1333+
"STLtools::fillSignedDistance is only available in 3D");
1334+
1335+
#if (AMREX_SPACEDIM != 3)
1336+
amrex::ignore_unused(this, mf, nghost, geom);
1337+
#else
1338+
this->fill(mf, nghost, geom, 1._rt, -1._rt);
1339+
1340+
const auto plo = geom.ProbLoArray();
1341+
const auto dx = geom.CellSizeArray();
1342+
1343+
auto ixt = mf.ixType();
1344+
RealVect offset(AMREX_D_DECL(ixt.cellCentered(0) ? 0.5_rt : 0.0_rt,
1345+
ixt.cellCentered(1) ? 0.5_rt : 0.0_rt,
1346+
ixt.cellCentered(2) ? 0.5_rt : 0.0_rt));
1347+
1348+
auto const& ma = mf.arrays();
1349+
1350+
if (m_bvh_optimization) {
1351+
auto const* bvh_root = m_bvh_nodes.data();
1352+
ParallelFor(mf, nghost, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
1353+
{
1354+
XDim3 coords = {plo[0]+(static_cast<Real>(i)+offset[0])*dx[0],
1355+
plo[1]+(static_cast<Real>(j)+offset[1])*dx[1],
1356+
plo[2]+(static_cast<Real>(k)+offset[2])*dx[2]};
1357+
auto d2 = bvh_d2(coords, bvh_root);
1358+
ma[b](i,j,k) *= std::sqrt(d2);
1359+
});
1360+
} else {
1361+
auto const* tri_pts = m_tri_pts_d.data();
1362+
int num_triangles = m_num_tri;
1363+
ParallelFor(mf, nghost, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
1364+
{
1365+
XDim3 coords = {plo[0]+(static_cast<Real>(i)+offset[0])*dx[0],
1366+
plo[1]+(static_cast<Real>(j)+offset[1])*dx[1],
1367+
plo[2]+(static_cast<Real>(k)+offset[2])*dx[2]};
1368+
auto d2 = std::numeric_limits<Real>::max();
1369+
for (int tr = 0; tr < num_triangles; ++tr) {
1370+
auto tmp = pt_tri_min_d2(coords, tri_pts[tr]);
1371+
d2 = std::min(d2, tmp);
1372+
}
1373+
ma[b](i,j,k) *= std::sqrt(d2);
1374+
});
1375+
}
1376+
Gpu::streamSynchronize();
1377+
#endif
1378+
}
1379+
11161380
}

0 commit comments

Comments
 (0)