Skip to content

Commit 1a61cd9

Browse files
weymouthweymouth
committed
Clean up closest() and add stl tests
Cleaned up a few bugs/confusing bits in closest(). This is not clever (!) but I think it is now correct. It only prunes when the BV (square) distance is greater than the current best distance. This is even applied at the leaf level since that distance function is much much faster than the triangle's. If you want a signed distance function, you need to initialize to floatmax (or the size of the total bounding box), forcing you to find the closest triangle so the sign is correct. However, when the mesh defines a membrane body (or you just want the unsigned distance) you can initialize with the max BDIM value. I've set this cutoff to 4 everywhere I noticed. A more efficient traversal will require a stack. I'm not positive about the best way to implement a stack on each thread of a GPU... and I'm also not sure how much better it will work. So let's try this out first.
1 parent b57f69e commit 1a61cd9

File tree

4 files changed

+5877
-45
lines changed

4 files changed

+5877
-45
lines changed

‎ext/WaterLilyGeometryBasicsExt.jl‎

Lines changed: 23 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ Adapt.@adapt_structure Meshbody
2323

2424
"""
2525
MeshBody(mesh::Union{Mesh, String};
26-
map::Function=(x,t)->x, boundary::Bool=true, half_thk::T=0.f0,
26+
map::Function=(x,t)->x, boundary::Bool=false, half_thk::T=0.f0,
2727
scale::T=1.f0, mem=Array, primitives::Union{BBox, BSphere}) where T
2828
2929
Constructor for a MeshBody:
@@ -57,7 +57,7 @@ using LinearAlgebra: cross
5757
import ImplicitBVH: BoundingVolume,BBox,BSphere
5858
@fastmath @inline inside(x::SVector, b::BoundingVolume) = inside(x, b.volume)
5959
@fastmath @inline inside(x::SVector, b::BBox) = all(b.lo.-4 .≤ x) && all(x .≤ b.up.+4)
60-
@fastmath @inline inside(x::SVector, b::BSphere) = sum(abs2,x .- b.x) - b.r^2 ≤ 4
60+
@fastmath @inline inside(x::SVector, b::BSphere) = √sum(abs2,x .- b.x) - b.r ≤ 4
6161

6262
import WaterLily: ×
6363
# linear shape function to interpolate inside element
@@ -66,25 +66,34 @@ import WaterLily: ×
6666
√sum(abs2,×(t[:,2]-t[:,1],p-t[:,1]))]
6767
@fastmath @inline get_velocity(p::SVector, tri, vel)= (dA=shape_value(p,tri); vel*dA/sum(dA))
6868

69+
# brute-force fallback when the BVH is not available
70+
closest(x::SVector,mesh) = findmin(tri->d²_fast(x, tri),mesh)
71+
6972
# traverse the BVH
7073
import ImplicitBVH: memory_index,unsafe_isvirtual
71-
@inline function closest(x::SVector{D,T},bvh::ImplicitBVH.BVH,mesh,::Val{true}) where {D,T}
74+
@inline function closest(x::SVector{D,T},bvh::ImplicitBVH.BVH,mesh;a=floatmax(T),verbose=false) where {D,T}
7275
tree = bvh.tree; length_nodes = length(bvh.nodes)
73-
u=Int32(0);a=d=T(64) # initial guess #TODO sensitive to initial a
76+
u=ncheck=lcheck=tcheck=Int32(0) # initialize
7477
# Depth-First-Search
75-
i=2; for _ in 1:4tree.levels^2 # prevent infinite loops
78+
i=1; while true
7679
@inbounds j = memory_index(tree,i)
7780
if j ≤ length_nodes # we are on a node
78-
inside(x, bvh.nodes[j]) && (i = 2i; continue) # go deeper
81+
verbose && (ncheck += 1)
82+
dist(x, bvh.nodes[j]) < a && (i = 2i; continue) # go deeper if closer than current best
7983
else # we reached a leaf
80-
@inbounds j = bvh.leaves[j-length_nodes].index # correct index in mesh
81-
d = d²_fast(x, mesh[j])
82-
d<a && (a=d; u=Int32(j)) # Replace current best
84+
verbose && (lcheck += 1)
85+
if dist(x, bvh.leaves[j-length_nodes]) < a
86+
verbose && (tcheck += 1)
87+
@inbounds j = bvh.leaves[j-length_nodes].index # correct index in mesh
88+
d = d²_fast(x, mesh[j])
89+
d<a && (a=d; u=Int32(j)) # Replace current best
90+
end
8391
end
8492
i = i>>trailing_ones(i)+1 # go to sibling, or uncle etc.
8593
(i==1 || unsafe_isvirtual(tree, i)) && break # search complete!
8694
end
87-
return u,a
95+
verbose && println("Checked $ncheck nodes, $lcheck leaves, $tcheck triangles")
96+
return a,u
8897
end
8998

9099
@inline function down(x,l,r,i)
@@ -94,28 +103,8 @@ end
94103
((dr ≤ 0) || (dr < dl)) && return 2i+1
95104
return 2i
96105
end
97-
@inline function closest(x::SVector{D,T},bvh::BVH,mesh,::Val{false}) where {D,T}
98-
tree = bvh.tree; length_nodes = length(bvh.nodes)
99-
u=Int32(1); leaf=0; a=d=T(1e6) # initial guess
100-
# start at top
101-
i=1; while leaf < 2 # max check two leafs
102-
@inbounds j = memory_index(tree, i)
103-
if j ≤ length_nodes # we are on a node
104-
@inbounds Il = memory_index(bvh.tree, 2i)
105-
@inbounds Ir = memory_index(bvh.tree, 2i+1)
106-
@inbounds i = down(x, bvh.nodes[Il], bvh.nodes[Ir], i)
107-
else # we reached a leaf
108-
@inbounds j = bvh.leaves[j-length_nodes].index # correct index in mesh
109-
d = d²_fast(x, mesh[j])
110-
d<a && (a=d; u=Int32(j)) # Replace current best
111-
i = i%2==0 ? i+1 : i-1; leaf += 1 # check sibling
112-
end
113-
unsafe_isvirtual(tree, i) && (i = i%2==0 ? i+1 : i-1) # go down sibling if virtual
114-
end
115-
return u,a
116-
end
117106

118-
# compute the distance to primitive
107+
# compute the square distance to primitive
119108
dist(x, b::BSphere) = sum(abs2,x .- b.x) - b.r
120109
function dist(x, b::BBox)
121110
c = (b.up .+ b.lo) ./ 2
@@ -124,16 +113,6 @@ function dist(x, b::BBox)
124113
end
125114
dist(x, b::BoundingVolume) = dist(x, b.volume)
126115

127-
# old brute-force function
128-
@inline function closest(x::SVector,mesh)
129-
u=Int32(1); a=b=d²_fast(x, mesh[1]) # fast method
130-
for I in 2:length(mesh)
131-
b = d²_fast(x, mesh[I])
132-
b<a && (a=b; u=I) # Replace current best
133-
end
134-
return u,a
135-
end
136-
137116
# locate the closest point p to x on triangle tri
138117
function locate(x::SVector{T},tri::SMatrix{T}) where T
139118
# unpack the triangle vertices
@@ -174,11 +153,10 @@ function WaterLily.measure(body::Meshbody,x::SVector{D,T},t;fastd²=Inf) where {
174153
# map to correct location
175154
ξ = body.map(x,t)
176155
# before we try the bvh
177-
!inside(ξ,body.bvh.nodes[1]) && return (T(8),zero(x),zero(x))
156+
!inside(ξ,body.bvh.nodes[1]) && return (T(4),zero(x),zero(x))
178157
# locate the point on the mesh
179-
u,d⁰ = closest(ξ,body.bvh,body.mesh,Val(true))
180-
u==0 && return (T(8),zero(x),zero(x)) # no closest found
181-
# u==0 && return check_inside(ξ, body.bvh) # no closest found
158+
_,u = closest(ξ,body.bvh,body.mesh;a = body.boundary ? floatmax(T) : T(16))
159+
u==0 && return (T(4),zero(x),zero(x)) # no triangles within distance "a"
182160
# compute the normal and distance
183161
n,p = normal(body.mesh[u]),SVector(locate(ξ,body.mesh[u]))
184162
# signed Euclidian distance

‎test/maintests.jl‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -750,4 +750,30 @@ end
750750
body += AutoBody((x,t)->42.f0) # the answer!
751751
@test GPUArrays.@allowscalar all(measure(body, x1.+1, 0) .≈ (0,[0,0,1],[1,1,1]))
752752
end
753+
754+
# Analytic SDF for a box with half-extents (a,b,c) centered at origin
755+
function box_sdf(x::SVector{3,T}, half_extents::SVector{3,T}=SA{T}[1, 1, 1]) where T
756+
q = abs.(x) .- half_extents
757+
return √sum(abs2,max.(q, zero(T))) + min(maximum(q), zero(T))
758+
end
759+
760+
# Load the box mesh using the public API
761+
box_file = joinpath(@__DIR__, "meshes", "box.stl")
762+
sphere_file = joinpath(@__DIR__, "meshes", "sphere.stl")
763+
for mem in arrays
764+
body = MeshBody(box_file; scale=1.f0, mem, boundary=true)
765+
for x in [SA{T}[0, 0, 0], SA{T}[0.5, 0.5, 0.5],
766+
SA{T}[1, 0, 0], SA{T}[0, 1, 0], SA{T}[0, 0, 1],
767+
SA{T}[2, 0, 0], SA{T}[1.5, 1.5, 1.5], SA{T}[1, 1, 1]]
768+
@show x, sdf(body, x, 0), box_sdf(x)
769+
@test sdf(body, x, 0) ≈ box_sdf(x)
770+
end
771+
body = MeshBody(sphere_file; scale=1.f0, mem, boundary=true)
772+
for x in [SA{T}[0, 0, 0], SA{T}[0.5, 0.5, 0.5],
773+
SA{T}[1, 0, 0], SA{T}[0, 1, 0], SA{T}[0, 0, 1],
774+
SA{T}[2, 0, 0], SA{T}[1.5, 1.5, 1.5], SA{T}[1, 1, 1]]
775+
@show x, sdf(body, x, 0), √sum(abs2, x)-0.5f0
776+
@test isapprox(sdf(body, x, 0), √sum(abs2, x)-0.5f0, rtol=15e-3)
777+
end
778+
end
753779
end

‎test/meshes/box.stl‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
solid vcg
2+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
3+
outer loop
4+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
5+
vertex 1.000000e+00 1.000000e+00 -1.000000e+00
6+
vertex 1.000000e+00 -1.000000e+00 -1.000000e+00
7+
endloop
8+
endfacet
9+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
10+
outer loop
11+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
12+
vertex -1.000000e+00 1.000000e+00 -1.000000e+00
13+
vertex 1.000000e+00 1.000000e+00 -1.000000e+00
14+
endloop
15+
endfacet
16+
facet normal -4.082483e-01 -4.082483e-01 8.164966e-01
17+
outer loop
18+
vertex -1.000000e+00 -1.000000e+00 1.000000e+00
19+
vertex 1.000000e+00 -1.000000e+00 1.000000e+00
20+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
21+
endloop
22+
endfacet
23+
facet normal -4.082483e-01 -4.082483e-01 8.164966e-01
24+
outer loop
25+
vertex -1.000000e+00 -1.000000e+00 1.000000e+00
26+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
27+
vertex -1.000000e+00 1.000000e+00 1.000000e+00
28+
endloop
29+
endfacet
30+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
31+
outer loop
32+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
33+
vertex 1.000000e+00 -1.000000e+00 -1.000000e+00
34+
vertex 1.000000e+00 -1.000000e+00 1.000000e+00
35+
endloop
36+
endfacet
37+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
38+
outer loop
39+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
40+
vertex 1.000000e+00 -1.000000e+00 1.000000e+00
41+
vertex -1.000000e+00 -1.000000e+00 1.000000e+00
42+
endloop
43+
endfacet
44+
facet normal -4.082483e-01 8.164966e-01 -4.082483e-01
45+
outer loop
46+
vertex -1.000000e+00 1.000000e+00 -1.000000e+00
47+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
48+
vertex 1.000000e+00 1.000000e+00 -1.000000e+00
49+
endloop
50+
endfacet
51+
facet normal -4.082483e-01 8.164966e-01 -4.082483e-01
52+
outer loop
53+
vertex -1.000000e+00 1.000000e+00 -1.000000e+00
54+
vertex -1.000000e+00 1.000000e+00 1.000000e+00
55+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
56+
endloop
57+
endfacet
58+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
59+
outer loop
60+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
61+
vertex -1.000000e+00 1.000000e+00 1.000000e+00
62+
vertex -1.000000e+00 1.000000e+00 -1.000000e+00
63+
endloop
64+
endfacet
65+
facet normal -5.773503e-01 -5.773503e-01 -5.773503e-01
66+
outer loop
67+
vertex -1.000000e+00 -1.000000e+00 -1.000000e+00
68+
vertex -1.000000e+00 -1.000000e+00 1.000000e+00
69+
vertex -1.000000e+00 1.000000e+00 1.000000e+00
70+
endloop
71+
endfacet
72+
facet normal 8.164966e-01 -4.082483e-01 -4.082483e-01
73+
outer loop
74+
vertex 1.000000e+00 -1.000000e+00 -1.000000e+00
75+
vertex 1.000000e+00 1.000000e+00 -1.000000e+00
76+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
77+
endloop
78+
endfacet
79+
facet normal 8.164966e-01 -4.082483e-01 -4.082483e-01
80+
outer loop
81+
vertex 1.000000e+00 -1.000000e+00 -1.000000e+00
82+
vertex 1.000000e+00 1.000000e+00 1.000000e+00
83+
vertex 1.000000e+00 -1.000000e+00 1.000000e+00
84+
endloop
85+
endfacet
86+
endsolid vcg

0 commit comments

Comments
 (0)