Add a nl_min_dist_sqr function

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

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)

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