diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index 2708e58..71446dd 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -131,6 +131,7 @@ export NeighborListInputs export NeighborLists export NeighborListEdge export is_neighbor_list_subset +export nl_min_dist_sqr export setup_neighbors_naive! export setup_neighbors_sort_sweep! @@ -145,6 +146,7 @@ let NL = NeighborLists{DefaultCollidePolicy, DefaultPairParams}, ) precompile(setup_neighbors_sort_sweep!, (NL, Pos, Inputs)) precompile(setup_neighbors_naive!, (NL, Pos, Inputs)) + precompile(nl_min_dist_sqr, (Pos, NL)) end end end diff --git a/src/neighbor-lists.jl b/src/neighbor-lists.jl index fceab96..f2e4b97 100644 --- a/src/neighbor-lists.jl +++ b/src/neighbor-lists.jl @@ -423,6 +423,56 @@ function is_neighbor_list_subset(subset::NeighborLists, superset::NeighborLists) return true end + +# Branch-free check for any NaN vertex coordinate in a simplex. +@inline function _simplex_has_nan(simplex) + bad = false + for v in simplex, x in v + bad |= isnan(x) + end + bad +end + +# Return the minimum squared distance between the object pairs in the neighbor +# list `nl`, `NaN` if any position used is `NaN`, or `Inf` if `nl` is empty. +function _nl_min_dist_sqr(pos, nl::AbstractVector{<:NeighborListEdge}) + T = eltype(eltype(pos)) + min_d2 = typemax(T) + # dist_sqr and min_fast have undefined NaN behavior, so track NaN separately + # by checking the input positions. + bad = false + for edge in nl + local a = load_positions(pos, edge.a) + local b = load_positions(pos, edge.b) + bad |= _simplex_has_nan(a) | _simplex_has_nan(b) + min_d2 = Base.FastMath.min_fast(min_d2, dist_sqr(a, b)) + end + ifelse(bad, T(NaN), min_d2) +end + +""" + nl_min_dist_sqr(pos, s::NeighborLists)::NamedTuple + +Return the minimum squared distance in each neighbor list in `s`, +keyed by the field names of `NeighborLists`. Empty lists have a minimum of `Inf`. +If any position used by a list is `NaN`, that list's minimum is `NaN`. + +Useful for checking if any objects are dangerously close, for example close +enough to pass through each other in one time step. +""" +function nl_min_dist_sqr(pos, s::NeighborLists) + (; + PPNL = _nl_min_dist_sqr(pos, s.PPNL), + PCNL = _nl_min_dist_sqr(pos, s.PCNL), + PLNL = _nl_min_dist_sqr(pos, s.PLNL), + PTNL = _nl_min_dist_sqr(pos, s.PTNL), + CCNL = _nl_min_dist_sqr(pos, s.CCNL), + CLNL = _nl_min_dist_sqr(pos, s.CLNL), + LLNL = _nl_min_dist_sqr(pos, s.LLNL), + ) +end + + function _prepare_neighbor_lists!(s::NeighborLists, inputs::NeighborListInputs) s.policy = inputs.policy empty!(s.PPNL) diff --git a/test/test-neighbor-lists.jl b/test/test-neighbor-lists.jl index cb93b2b..53fedc8 100644 --- a/test/test-neighbor-lists.jl +++ b/test/test-neighbor-lists.jl @@ -320,5 +320,44 @@ using Random end end end + @testset "nl_min_dist_sqr" begin + s = NeighborLists(policy) + pos = [ + SVector{3, Float32}(0, 0, 0), + SVector{3, Float32}(1, 0, 0), + SVector{3, Float32}(3, 0, 0), + SVector{3, Float32}(0, 0, 2), + SVector{3, Float32}(1, 0, 2), + ] + params = DefaultPairParams(1.0f0) + # Empty lists have an Inf minimum + @test nl_min_dist_sqr(pos, s) === (; + PPNL=Inf32, PCNL=Inf32, PLNL=Inf32, PTNL=Inf32, + CCNL=Inf32, CLNL=Inf32, LLNL=Inf32, + ) + push!(s.PPNL, NeighborListEdge(PointIdxPart(1), PointIdxPart(3), 1.0f0, params)) + push!(s.PPNL, NeighborListEdge(PointIdxPart(1), PointIdxPart(2), 1.0f0, params)) + push!(s.LLNL, NeighborListEdge(LineIdxPart(1, 2), LineIdxPart(4, 5), 1.0f0, params)) + min_d2s = nl_min_dist_sqr(pos, s) + @test min_d2s.PPNL === 1.0f0 + @test min_d2s.LLNL === 4.0f0 + @test min_d2s.PTNL === Inf32 + # Works with reinterpreted Float64 positions + pos64 = reinterpret(SVector{3, Float64}, vec(Float64.(stack(pos)))) + min_d2s64 = nl_min_dist_sqr(pos64, s) + @test min_d2s64.PPNL === 1.0 + @test min_d2s64.LLNL === 4.0 + @test min_d2s64.PCNL === Inf + # NaN positions result in NaN, regardless of position in the list + for nan_idx in (1, 2, 3) + nan_pos = copy(pos) + nan_pos[nan_idx] = SVector{3, Float32}(NaN32, 0, 0) + @test nl_min_dist_sqr(nan_pos, s).PPNL === NaN32 + end + nan_pos = copy(pos) + nan_pos[4] = SVector{3, Float32}(NaN32, 0, 0) + @test nl_min_dist_sqr(nan_pos, s).PPNL === 1.0f0 + @test nl_min_dist_sqr(nan_pos, s).LLNL === NaN32 + end end nothing