From 78ccfafe7a344e17ad45178478808d90d189631e Mon Sep 17 00:00:00 2001 From: nhz2 Date: Tue, 11 Aug 2026 09:51:44 +0200 Subject: [PATCH 1/4] Add collide forces --- Project.toml | 10 +- README.md | 47 ++++- src/SimplexCellLists.jl | 18 ++ src/collide-forces.jl | 353 ++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 + test/test-collide-forces.jl | 177 ++++++++++++++++++ 6 files changed, 602 insertions(+), 5 deletions(-) create mode 100644 src/collide-forces.jl create mode 100644 test/test-collide-forces.jl diff --git a/Project.toml b/Project.toml index 1156911..1880fcf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,18 +1,20 @@ name = "SimplexCellLists" uuid = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4" -authors = ["Nathan Zimmerberg"] version = "0.2.1-dev" +authors = ["Nathan Zimmerberg"] + +[workspace] +projects = ["test", "benchmark"] [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" +ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "2" +ChunkSplitters = "3" LinearAlgebra = "1" StaticArrays = "1" julia = "1.12" - -[workspace] -projects = ["test", "benchmark"] diff --git a/README.md b/README.md index 0ec8e24..bc19b83 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ This Julia package contains data structures and algorithms for doing computations on pairs of 3D points, line segments, and triangles within a cutoff distance. -It provides cell lists for points ([`PointCellList`](src/pointcelllist.jl)) and line segments ([`LineSegCellList`](src/linesegcelllist.jl)) with fast nearby-neighbor mapping, and squared-distance functions (`dist_sqr`) for all pairs of points, line segments, and triangles. +It provides cell lists for points ([`PointCellList`](src/pointcelllist.jl)) and line segments ([`LineSegCellList`](src/linesegcelllist.jl)) with fast nearby-neighbor mapping, squared-distance functions (`dist_sqr`) for all pairs of points, line segments, and triangles, neighbor list construction ([`neighbor-lists.jl`](src/neighbor-lists.jl)) with a sort and sweep broad phase, and soft-sphere collide forces and energy ([`collide-forces.jl`](src/collide-forces.jl)) computed from the neighbor lists. This package is largely inspired by [CellListMap.jl](https://github.com/m3g/CellListMap.jl). @@ -29,3 +29,48 @@ end ``` `LineSegCellList`, `cell_line_seg_add!`, and `map_nearby_line_segs` are the line segment analogs, and `cell_points_clear!`/`cell_line_segs_clear!` reset a list for reuse. + +## Collide forces + +Neighbor lists and collide forces are configured with a `CollisionPolicy`. +The policy decides which objects and pairs participate in the neighbor lists +(`filter_object`, `filter_pair`), how per-object parameters combine into +per-pair parameters (`mix_params`), and which force law applies to each edge +(`nl_edge_forces!`). + +`DefaultCollisionPolicy` stores a stiffness and collision layer masks per +object, and applies a soft repulsive potential `E = k/2 * (L - d)²` when two +objects overlap, where `L` is the sum of the two radii and `d` is the closest +distance between them: + +```julia +using SimplexCellLists, StaticArrays + +policy = DefaultCollisionPolicy() + +# Two overlapping spheres: radius 0.5, centers 0.5 apart +pos = [SA[0.0, 0.0, 0.0], SA[0.5, 0.0, 0.0]] +params = DefaultObjectParams(10.0f0, UInt32(1), UInt32(0)) # stiffness, layers, no collide mask +inputs = NeighborListInputs(policy; + points = [PointIdxPart(1), PointIdxPart(2)], + p_radius = [0.5f0, 0.5f0], + p_params = [params, params], +) +nl = NeighborLists(policy) +setup_neighbors_sort_sweep!(nl, pos, inputs) + +force_energy = ForceEnergyFloat64(length(pos)) +collide_forces!(force_energy, pos, nl, Float64) +get_energy(force_energy) # 0.625 +get_force(force_energy) # [-2.5, 0.0, 0.0], [2.5, 0.0, 0.0] +``` + +`NeighborListInputs` also accepts line segments (`clines`/`lines`) and +triangles, `no_collide_pairs` exclusions, and a `skin` distance so lists can be +reused across steps. `setup_neighbors_naive!` is a reference implementation of +`setup_neighbors_sort_sweep!` useful for testing. + +A custom policy subtypes `CollisionPolicy{ObjectParams, PairParams}` with its +own parameter types and methods for `filter_object`, `filter_pair`, and +`mix_params`, and can define its own `nl_edge_forces!` methods to change the +force law. diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index 053a369..9fe4139 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -3,6 +3,7 @@ module SimplexCellLists using LinearAlgebra: LinearAlgebra, cross, dot, ⋅ using StaticArrays: StaticArrays, @SVector, SA, SVector using ArgCheck: @argcheck +using ChunkSplitters: chunks const Vec3 = SVector{3} const Simplex{N,T} = SVector{N, Vec3{T}} @@ -127,4 +128,21 @@ let NL = NeighborLists{DefaultCollisionPolicy, DefaultPairParams}, end end +include("collide-forces.jl") +export collide_forces! + +# Precompile the collide forces for the default policy +let NL = NeighborLists{DefaultCollisionPolicy, DefaultPairParams} + for FE in (ForceEnergyFloat64, ForceEnergyFixedPoint{30, 30}, ForceNoEnergyFixedPoint{30}) + for T in (Float32, Float64) + for Pos in ( + Vector{SVector{3, T}}, + typeof(reinterpret(SVector{3, T}, T[])), + ) + precompile(collide_forces!, (FE, Pos, NL, Type{Float64})) + end + end + end +end + end diff --git a/src/collide-forces.jl b/src/collide-forces.jl new file mode 100644 index 0000000..5963031 --- /dev/null +++ b/src/collide-forces.jl @@ -0,0 +1,353 @@ +""" + nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge, policy::CollisionPolicy, calc_type)::Nothing + +Add the collision force and energy of a single neighbor-list edge to `force_energy`. + +Methods dispatch on the index part types of `edge` and on `policy`, so a custom +`CollisionPolicy` with its own `PairParams` can define its own force laws. +""" +function nl_edge_forces! end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, PointIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + x1 = pos[edge.a.i] + x2 = pos[edge.b.i] + d = map(calc_type, x1 - x2) + d2cl = sum(abs2, d) + L0 = calc_type(edge.L) + if d2cl < L0^2 + L = _sqrt_fast(d2cl) + k = calc_type(edge.params.k) + ΔL = L0 - L + E = k*calc_type(1//2)*ΔL^2 + F = (k*ΔL*inv(L))*d + add_bead_force!(force_energy, edge.a.i, F) + add_bead_force!(force_energy, edge.b.i, -F) + add_energy!(force_energy, E) + end + nothing +end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + x1 = pos[edge.a.i] + y1 = pos[edge.b.i] + y2 = pos[edge.b.i + UInt32(1)] + r = map(calc_type, y2 - y1) + p = map(calc_type, x1 - y1) + t = clamp01nan((p⋅r)*inv(r⋅r)) + d = t*r - p + d2cl = sum(abs2, d) + L0 = calc_type(edge.L) + if d2cl < L0^2 + L = _sqrt_fast(d2cl) + k = calc_type(edge.params.k) + ΔL = L0 - L + E = k*calc_type(1//2)*ΔL^2 + F = (k*ΔL*inv(L))*d + add_bead_force!(force_energy, edge.a.i, -F) + add_bead_force!(force_energy, edge.b.i, (oneunit(t)-t)*F) + add_bead_force!(force_energy, edge.b.i + UInt32(1), t*F) + add_energy!(force_energy, E) + end + nothing +end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + x1 = pos[edge.a.i] + y1 = pos[edge.b.i] + y2 = pos[edge.b.j] + r = map(calc_type, y2 - y1) + p = map(calc_type, x1 - y1) + t = clamp01nan((p⋅r)*inv(r⋅r)) + d = t*r - p + d2cl = sum(abs2, d) + L0 = calc_type(edge.L) + if d2cl < L0^2 + L = _sqrt_fast(d2cl) + k = calc_type(edge.params.k) + ΔL = L0 - L + E = k*calc_type(1//2)*ΔL^2 + F = (k*ΔL*inv(L))*d + add_bead_force!(force_energy, edge.a.i, -F) + add_bead_force!(force_energy, edge.b.i, (oneunit(t)-t)*F) + add_bead_force!(force_energy, edge.b.j, t*F) + add_energy!(force_energy, E) + end + nothing +end + +""" +Core helper function for line-line collide force calculation. +Returns (do_force::Bool, u, fp1, fq0, fq1) where: +- do_force: whether forces should be applied +- u: energy (without k factor) +- fp1: force on P1 (without k factor) +- fq0: force on Q0 (without k factor) +- fq1: force on Q1 (without k factor) +Force on P0 can be computed as -(fp1+fq0+fq1) +""" +@inline function line_line_force_core(P0, P1, Q0, Q1, calc_type, switchover_scale_unitless, L0) + P = map(calc_type, P1-P0) + Q = map(calc_type, Q1-Q0) + P0mQ0 = map(calc_type, P0-Q0) + a = P ⋅ P + b = P ⋅ Q + c = Q ⋅ Q + d = P ⋅ P0mQ0 + e = Q ⋅ P0mQ0 + f = P0mQ0 ⋅ P0mQ0 + Δ = a*c - b^2 + #assuming both segments are not zero length + #critical points + s0 = clamp01nan(-d/a) + s1 = clamp01nan((b-d)/a) + t0 = clamp01nan(e/c) + t1 = clamp01nan((b+e)/c) + sbar = clamp01nan((b*e - c*d)/Δ) + tbar = clamp01nan((a*e - b*d)/Δ) + r(s,t) = a*s^2 - 2b*s*t + c*t^2 + 2d*s - 2e*t + mins = sbar + mint = tbar + minr = r(sbar,tbar) + if r(s0,0) < minr + minr = r(s0,0) + mins = s0 + mint = zero(a) + end + if r(s1,1) < minr + minr = r(s1,1) + mins = s1 + mint = one(a) + end + if r(0,t0) < minr + minr = r(0,t0) + mins = zero(a) + mint = t0 + end + if r(1,t1) < minr + minr = r(1,t1) + mins = one(a) + mint = t1 + end + d2cl = relu(f+minr) + if d2cl < L0^2 + #return the energy and force on each point give s and t + # doesn't have the e.k factor + function u_f(s,t) + local Pcl = P0mQ0 + s*P + local Qcl = t*Q + local Rcl = Qcl - Pcl + local L = norm_fast(Rcl) + local ΔL = relu(L0 - L) + local u = calc_type(1//2)*ΔL^2 + local FP = (-ΔL*inv(L))*Rcl + (u, s*FP, -(1-t)*FP, -t*FP) + end + switchover_scale = calc_type(switchover_scale_unitless)*inv(L0)^2 # units of 1/nm^2, + switchover = tanh(Δ/sqrt(a*c)*switchover_scale) + #P0 is set as origin for now, subtract forces later to get final fp0 + ucl, fp1cl, fq0cl, fq1cl = u_f(mins,mint) + us0, fp1s0, fq0s0, fq1s0 = u_f(s0,0) + us1, fp1s1, fq0s1, fq1s1 = u_f(s1,1) + ut0, fp1t0, fq0t0, fq1t0 = u_f(0,t0) + ut1, fp1t1, fq0t1, fq1t1 = u_f(1,t1) + uother = +( + us0, + us1, + ut0, + ut1, + ) + fp1other = +( + fp1s0, + fp1s1, + fp1t0, + fp1t1, + ) + fq0other = +( + fq0s0, + fq0s1, + fq0t0, + fq0t1, + ) + fq1other = +( + fq1s0, + fq1s1, + fq1t0, + fq1t1, + ) + u = switchover*ucl + calc_type(1//2)*(1-switchover)*uother + p = sqrt(a) + q = sqrt(c) + g = b/(p*q) + uswitched = (ucl + -calc_type(1//2)*uother) + fp1 = switchover*fp1cl + calc_type(1//2)*(1-switchover)*fp1other + d_switchover_dp1 = -switchover_scale*(1-switchover^2)*((q/p+b*g/a)*P - 2*g*Q) + fp1 += d_switchover_dp1*uswitched + + fq0 = switchover*fq0cl + calc_type(1//2)*(1-switchover)*fq0other + fq1 = switchover*fq1cl + calc_type(1//2)*(1-switchover)*fq1other + d_switchover_dq0 = switchover_scale*(1-switchover^2)*((p/q+b*g/c)*Q - 2*g*P) + d_switchover_dq1 = - d_switchover_dq0 + fq0 += d_switchover_dq0*uswitched + fq1 += d_switchover_dq1*uswitched + + return true, u, fp1, fq0, fq1 + else + return false, zero(L0), zero(P), zero(P), zero(P) + end +end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + P0 = pos[edge.a.i] + P1 = pos[edge.a.i + UInt32(1)] + Q0 = pos[edge.b.i] + Q1 = pos[edge.b.i + UInt32(1)] + L0 = calc_type(edge.L) + do_force, u, fp1, fq0, fq1 = line_line_force_core(P0, P1, Q0, Q1, calc_type, policy.switchover_scale_unitless, L0) + if do_force + ke = calc_type(edge.params.k) + add_bead_force!(force_energy, edge.a.i, -ke*(fp1+fq0+fq1)) + add_bead_force!(force_energy, edge.a.i + UInt32(1), ke*fp1) + + add_bead_force!(force_energy, edge.b.i, ke*fq0) + add_bead_force!(force_energy, edge.b.i + UInt32(1), ke*fq1) + + add_energy!(force_energy, u*ke) + end + nothing +end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + P0 = pos[edge.a.i] + P1 = pos[edge.a.i + UInt32(1)] + Q0 = pos[edge.b.i] + Q1 = pos[edge.b.j] + L0 = calc_type(edge.L) + do_force, u, fp1, fq0, fq1 = line_line_force_core(P0, P1, Q0, Q1, calc_type, policy.switchover_scale_unitless, L0) + if do_force + ke = calc_type(edge.params.k) + add_bead_force!(force_energy, edge.a.i, -ke*(fp1+fq0+fq1)) + add_bead_force!(force_energy, edge.a.i + UInt32(1), ke*fp1) + + add_bead_force!(force_energy, edge.b.i, ke*fq0) + add_bead_force!(force_energy, edge.b.j, ke*fq1) + + add_energy!(force_energy, u*ke) + end + nothing +end + +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{LineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + P0 = pos[edge.a.i] + P1 = pos[edge.a.j] + Q0 = pos[edge.b.i] + Q1 = pos[edge.b.j] + L0 = calc_type(edge.L) + do_force, u, fp1, fq0, fq1 = line_line_force_core(P0, P1, Q0, Q1, calc_type, policy.switchover_scale_unitless, L0) + if do_force + ke = calc_type(edge.params.k) + add_bead_force!(force_energy, edge.a.i, -ke*(fp1+fq0+fq1)) + add_bead_force!(force_energy, edge.a.j, ke*fp1) + + add_bead_force!(force_energy, edge.b.i, ke*fq0) + add_bead_force!(force_energy, edge.b.j, ke*fq1) + + add_energy!(force_energy, u*ke) + end + nothing +end + +#= +Based on https://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf +// David Eberly, Geometric Tools, Redmond WA 98052 +// Copyright (c) 1998-2022 +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt +// Version: 6.0.2022.01.06 +=# +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, TriangleIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) + P = pos[edge.a.i] + B = pos[edge.b.i] + E0 = map(calc_type, pos[edge.b.j] - B) + E1 = map(calc_type, pos[edge.b.k] - B) + BP = map(calc_type, B - P) + a = E0 ⋅ E0 + b = E0 ⋅ E1 + c = E1 ⋅ E1 + d = E0 ⋅ BP + e = E1 ⋅ BP + f = BP ⋅ BP + Δ = a*c - b^2 + invΔ = inv(Δ) + sbar = clamp01nan((b*e - c*d)*invΔ) + tbar = clamp01nan((-a*e + b*d)*invΔ) + + # this ensures sbar and tbar are in the domain + # if they are out of the domain, the min distance is on a boundary + outside = sbar+tbar > one(sbar) + sbar = ifelse(outside,zero(sbar),sbar) + tbar = ifelse(outside,zero(tbar),tbar) + s0 = clamp01nan(-d/a) + t0 = clamp01nan(-e/c) + sd = clamp01nan(-(b-c+d-e)/(a-2b+c)) + td = one(sd) - sd + r(s,t) = a*s^2 + 2b*s*t + c*t^2 + 2d*s + 2e*t + + mins = sbar + mint = tbar + minr = r(sbar,tbar) + if r(s0,0) < minr + minr = r(s0,0) + mins = s0 + mint = zero(a) + end + if r(0,t0) < minr + minr = r(0,t0) + mins = zero(a) + mint = t0 + end + if r(sd,td) < minr + minr = r(sd,td) + mins = sd + mint = td + end + d2cl = relu(f+minr) + L0 = calc_type(edge.L) + if d2cl < L0^2 + # dvec is the vector from P to closest point on triangle + dvec = BP + mins*E0 + mint*E1 + L = _sqrt_fast(d2cl) + k = calc_type(edge.params.k) + ΔL = L0 - L + E = k*calc_type(1//2)*ΔL^2 + F = (k*ΔL*inv(L))*dvec + add_bead_force!(force_energy, edge.a.i, -F) + # Forces on triangle vertices + # Closest point on triangle is B + mins*E0 + mint*E1 = (1-mins-mint)*B + mins*(B+E0) + mint*(B+E1) + # = (1-mins-mint)*y[1] + mins*y[2] + mint*y[3] + w0 = oneunit(mins) - mins - mint + add_bead_force!(force_energy, edge.b.i, w0*F) + add_bead_force!(force_energy, edge.b.j, mins*F) + add_bead_force!(force_energy, edge.b.k, mint*F) + add_energy!(force_energy, E) + end + nothing +end + +function nl_forces!(force_energy::ForceEnergy, pos, nl, policy::CollisionPolicy, calc_type::T, chunk, nthreads) where T + chunk > length(nl) && return + for edge in chunks(nl; n=nthreads)[chunk] + nl_edge_forces!(force_energy, pos, edge, policy, calc_type) + end +end + +function collide_forces!(force_energy::ForceEnergy, pos, s::NeighborLists, calc_type; chunk=1, nthreads=1) + nl_forces!(force_energy, pos, s.PPNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.PCNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.PLNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.PTNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.CCNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.CLNL, s.policy, calc_type, chunk, nthreads) + nl_forces!(force_energy, pos, s.LLNL, s.policy, calc_type, chunk, nthreads) + nothing +end diff --git a/test/runtests.jl b/test/runtests.jl index 6fff2b4..b7d413d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,3 +12,5 @@ include("test-force-energy.jl") include("test-neighbor-lists.jl") +include("test-collide-forces.jl") + diff --git a/test/test-collide-forces.jl b/test/test-collide-forces.jl new file mode 100644 index 0000000..36f5af7 --- /dev/null +++ b/test/test-collide-forces.jl @@ -0,0 +1,177 @@ +using Test +using StaticArrays +using Rotations: RotMatrix +using SimplexCellLists +using SimplexCellLists: dist_sqr, nl_edge_forces!, NeighborListEdge + +flat_force(force_energy) = reinterpret(Float64, get_force(force_energy)) + +""" +Error if the collide force and energy violates a property +`a` and `b` are vectors of three vectors +""" +function check_collide_props(a, b) + policy = DefaultCollisionPolicy() + positions = SVector{3, Float64}[] + a_i, b_i = map((a,b)) do y + local idxs = UInt32[] + if rand(Bool) + # Add padding + push!(positions, SA[NaN,NaN,NaN]) + end + for x in y + while rand(Bool) + # Add padding + push!(positions, SA[NaN,NaN,NaN]) + end + push!(positions, x) + push!(idxs, UInt32(length(positions))) + end + if rand(Bool) + # Add padding + push!(positions, SA[NaN,NaN,NaN]) + end + if length(y) == 1 + PointIdxPart(idxs...) + elseif length(idxs) == 2 + if rand(Bool) || idxs[2] != idxs[1] + 1 + LineIdxPart(idxs...) + else + CLineIdxPart(idxs[1]) + end + elseif length(idxs) == 3 + TriangleIdxPart(idxs...) + else + error("unsupported shape") + end + end + # flip order if needed + if a_i isa LineIdxPart && b_i isa CLineIdxPart + a_i, b_i = b_i, a_i + end + k = 10f0 + dcl = sqrt(dist_sqr(SVector((a...,)), SVector((b...,)))) + # Out of range should have zero energy + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, Float32(dcl*0.9), DefaultPairParams(k)), policy, Float64) + @test iszero(get_energy(force_energy)) + @test all(iszero, get_force(force_energy)) + + # in range should have non zero energy proportional to k + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, Float32(dcl*1.5), DefaultPairParams(k)), policy, Float64) + e1 = get_energy(force_energy) + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, Float32(dcl*1.5), DefaultPairParams(10*k)), policy, Float64) + e2 = get_energy(force_energy) + @test e1 > 0 + @test 10*e1 ≈ e2 + + # Finite difference test + L = Float32(dcl*1.5) + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + e0 = get_energy(force_energy) + f = flat_force(force_energy) + for i in 1:length(positions)*3 + local new_pos = copy(positions) + local Δx = 1E-5 + reinterpret(Float64, new_pos)[i] += Δx + local force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, new_pos, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + local Δe = get_energy(force_energy) - e0 + @test (f[i]*Δx + Δe) < 1.5E-7 + end + + # Translation invariance + L = Float32(dcl*1.5) + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + e0 = get_energy(force_energy) + f = flat_force(force_energy) + new_pos = positions .+ (SA[3.0, 2.0, 1.0],) + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, new_pos, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + @test e0 ≈ get_energy(force_energy) + @test f ≈ flat_force(force_energy) + + # Rotation invariance + L = Float32(dcl*1.5) + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + e0 = get_energy(force_energy) + f = flat_force(force_energy) + new_pos = (rand(RotMatrix{3}),) .* positions + force_energy = ForceEnergyFloat64(length(positions)) + nl_edge_forces!(force_energy, new_pos, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) + @test e0 ≈ get_energy(force_energy) + return +end + + +@testset "collide forces" begin + x = SA[1.0, 0.0, 0.0] + y = SA[0.0, 1.0, 0.0] + z = SA[0.0, 0.0, 1.0] + check_collide_props( + SA[0z], + SA[x], + ) + check_collide_props( + SA[z], + SA[x], + ) + check_collide_props( + SA[0z], + SA[x, 2x], + ) + check_collide_props( + SA[0z], + SA[x, 1.001x], + ) + check_collide_props( + SA[1.002x], + SA[x, 1.001x], + ) + check_collide_props( + SA[2y], + SA[-6x, 6x], + ) + for i in 1:5000 + check_collide_props( + SA[0.1*y + randn()*x + randn()*z], + SA[randn()*x + randn()*z, randn()*x + randn()*z, randn()*x + randn()*z], + ) + end + check_collide_props( + SA[0*x, -1x], + SA[x, 2x], + ) + check_collide_props( + SA[0y, -1y], + SA[x, 2x], + ) + check_collide_props( + SA[0.1y+x, 0.1y+2x], + SA[x, 2x], + ) + check_collide_props( + SA[0.1y+2x, 0.1y+x], + SA[x, 2x], + ) + check_collide_props( + SA[0.1y-2z, 0.1y+2z], + SA[-2x, 2x], + ) + check_collide_props( + SA[20y-2z, 20y+2z], + SA[-2x, 2x], + ) + for i in 1:5000 + check_collide_props( + SA[0.1*y + randn()*x + randn()*z, 0.1*y + randn()*x + randn()*z], + SA[randn()*x + randn()*z, randn()*x + randn()*z], + ) + end +end +nothing From a1afc4d4ec53dc6f0f53d297493cdcf03fe03567 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Tue, 11 Aug 2026 22:26:17 +0200 Subject: [PATCH 2/4] docs improvement --- README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index bc19b83..e5d322c 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,12 @@ This Julia package contains data structures and algorithms for doing computations on pairs of 3D points, line segments, and triangles within a cutoff distance. -It provides cell lists for points ([`PointCellList`](src/pointcelllist.jl)) and line segments ([`LineSegCellList`](src/linesegcelllist.jl)) with fast nearby-neighbor mapping, squared-distance functions (`dist_sqr`) for all pairs of points, line segments, and triangles, neighbor list construction ([`neighbor-lists.jl`](src/neighbor-lists.jl)) with a sort and sweep broad phase, and soft-sphere collide forces and energy ([`collide-forces.jl`](src/collide-forces.jl)) computed from the neighbor lists. +It provides: + +- Cell lists for points ([`PointCellList`](src/pointcelllist.jl)) and line segments ([`LineSegCellList`](src/linesegcelllist.jl)) with fast nearby-neighbor mapping. +- Squared-distance functions (`dist_sqr`) for all pairs of points, line segments, and triangles. +- Neighbor list construction ([`neighbor-lists.jl`](src/neighbor-lists.jl)) with a sort and sweep broad phase. +- Soft-sphere collide forces and energy ([`collide-forces.jl`](src/collide-forces.jl)) computed from the neighbor lists. This package is largely inspired by [CellListMap.jl](https://github.com/m3g/CellListMap.jl). From 37d491a319833744a611d248f903259948a83df9 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Tue, 11 Aug 2026 23:01:19 +0200 Subject: [PATCH 3/4] more tests --- README.md | 4 ++- src/collide-forces.jl | 19 ++++++++++-- test/test-collide-forces.jl | 61 +++++++++++++++++++++++++++++++++++-- 3 files changed, 78 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index e5d322c..ad0bb70 100644 --- a/README.md +++ b/README.md @@ -46,7 +46,9 @@ per-pair parameters (`mix_params`), and which force law applies to each edge `DefaultCollisionPolicy` stores a stiffness and collision layer masks per object, and applies a soft repulsive potential `E = k/2 * (L - d)²` when two objects overlap, where `L` is the sum of the two radii and `d` is the closest -distance between them: +distance between them. The pair stiffness `k` mixes the two per-object +stiffnesses like springs in series, `k = k₁*k₂/(k₁ + k₂)`, so the two +`10.0f0` stiffnesses in the example below mix to `k = 5`: ```julia using SimplexCellLists, StaticArrays diff --git a/src/collide-forces.jl b/src/collide-forces.jl index 5963031..6d137fd 100644 --- a/src/collide-forces.jl +++ b/src/collide-forces.jl @@ -1,7 +1,7 @@ """ nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge, policy::CollisionPolicy, calc_type)::Nothing -Add the collision force and energy of a single neighbor-list edge to `force_energy`. +Add the collide force and energy of a single neighbor-list edge to `force_energy`. Methods dispatch on the index part types of `edge` and on `policy`, so a custom `CollisionPolicy` with its own `PairParams` can define its own force laws. @@ -142,7 +142,7 @@ Force on P0 can be computed as -(fp1+fq0+fq1) local FP = (-ΔL*inv(L))*Rcl (u, s*FP, -(1-t)*FP, -t*FP) end - switchover_scale = calc_type(switchover_scale_unitless)*inv(L0)^2 # units of 1/nm^2, + switchover_scale = calc_type(switchover_scale_unitless)*inv(L0)^2 # units of 1/length^2 switchover = tanh(Δ/sqrt(a*c)*switchover_scale) #P0 is set as origin for now, subtract forces later to get final fp0 ucl, fp1cl, fq0cl, fq1cl = u_f(mins,mint) @@ -341,6 +341,21 @@ function nl_forces!(force_energy::ForceEnergy, pos, nl, policy::CollisionPolicy, end end +""" + collide_forces!(force_energy::ForceEnergy, pos, s::NeighborLists, calc_type; chunk=1, nthreads=1)::Nothing + +Add the collide forces and energy of every neighbor-list edge in `s` to +`force_energy`, doing the calculations in the floating-point type `calc_type`. + +Each edge's contribution is computed by [`nl_edge_forces!`](@ref), so a custom +`CollisionPolicy` can change the force laws. + +For multithreading, split the work into `nthreads` chunks: thread `t` calls +`collide_forces!` with `chunk=t` and the same `nthreads`, accumulating into +its own `force_energy`, and the per-thread accumulators are merged afterwards +with [`combine_force_energy!`](@ref). Every edge is visited exactly once +across the chunks. +""" function collide_forces!(force_energy::ForceEnergy, pos, s::NeighborLists, calc_type; chunk=1, nthreads=1) nl_forces!(force_energy, pos, s.PPNL, s.policy, calc_type, chunk, nthreads) nl_forces!(force_energy, pos, s.PCNL, s.policy, calc_type, chunk, nthreads) diff --git a/test/test-collide-forces.jl b/test/test-collide-forces.jl index 36f5af7..02a6132 100644 --- a/test/test-collide-forces.jl +++ b/test/test-collide-forces.jl @@ -80,7 +80,7 @@ function check_collide_props(a, b) local force_energy = ForceEnergyFloat64(length(positions)) nl_edge_forces!(force_energy, new_pos, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) local Δe = get_energy(force_energy) - e0 - @test (f[i]*Δx + Δe) < 1.5E-7 + @test abs(f[i]*Δx + Δe) < 1.5E-7 end # Translation invariance @@ -100,11 +100,13 @@ function check_collide_props(a, b) force_energy = ForceEnergyFloat64(length(positions)) nl_edge_forces!(force_energy, positions, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) e0 = get_energy(force_energy) - f = flat_force(force_energy) - new_pos = (rand(RotMatrix{3}),) .* positions + f0 = get_force(force_energy) + R = rand(RotMatrix{3}) + new_pos = (R,) .* positions force_energy = ForceEnergyFloat64(length(positions)) nl_edge_forces!(force_energy, new_pos, NeighborListEdge(a_i, b_i, L, DefaultPairParams(k)), policy, Float64) @test e0 ≈ get_energy(force_energy) + @test (R,) .* f0 ≈ get_force(force_energy) return end @@ -174,4 +176,57 @@ end ) end end + +@testset "collide_forces! integration" begin + policy = DefaultCollisionPolicy() + # 3×3×4 grid, spacing 0.4: with radius 0.3 overlaps are guaranteed, but + # no two objects touch exactly, so all forces stay finite + positions = [0.4*SA[mod(i-1, 3), mod((i-1)÷3, 3), (i-1)÷9] for i in 1:30] + params = DefaultObjectParams(10.0f0, UInt32(1), UInt32(0)) + inputs = NeighborListInputs(policy; + points = [PointIdxPart(i) for i in 1:10], + p_radius = fill(0.3f0, 10), + p_params = fill(params, 10), + clines = [CLineIdxPart(i) for i in 11:2:17], + c_radius = fill(0.3f0, 4), + c_params = fill(params, 4), + lines = [LineIdxPart(19, 20), LineIdxPart(21, 22)], + l_radius = fill(0.3f0, 2), + l_params = fill(params, 2), + triangles = [TriangleIdxPart(23, 24, 25)], + t_radius = [0.5f0], + t_params = [params], + ) + nl = NeighborLists(policy) + setup_neighbors_sort_sweep!(nl, positions, inputs) + # every edge type participates + for list in (nl.PPNL, nl.PCNL, nl.PLNL, nl.PTNL, nl.CCNL, nl.CLNL, nl.LLNL) + @test !isempty(list) + end + + fe = ForceEnergyFloat64(length(positions)) + collide_forces!(fe, positions, nl, Float64) + @test get_energy(fe) > 0 + @test all(f -> all(isfinite, f), get_force(fe)) + + # matches summing nl_edge_forces! over every edge + fe_manual = ForceEnergyFloat64(length(positions)) + for list in (nl.PPNL, nl.PCNL, nl.PLNL, nl.PTNL, nl.CCNL, nl.CLNL, nl.LLNL) + for edge in list + nl_edge_forces!(fe_manual, positions, edge, policy, Float64) + end + end + @test get_energy(fe_manual) == get_energy(fe) + @test get_force(fe_manual) == get_force(fe) + + # the chunked path visits every edge exactly once + for nthreads in (1, 2, 3, 7, 100) + local fe_chunked = ForceEnergyFloat64(length(positions)) + for chunk in 1:nthreads + collide_forces!(fe_chunked, positions, nl, Float64; chunk, nthreads) + end + @test get_energy(fe_chunked) ≈ get_energy(fe) + @test get_force(fe_chunked) ≈ get_force(fe) + end +end nothing From 787f4362c90d8adbdbe9605e8ba7c6c43681f0b7 Mon Sep 17 00:00:00 2001 From: nhz2 Date: Tue, 11 Aug 2026 23:43:56 +0200 Subject: [PATCH 4/4] fix allocations --- Project.toml | 2 -- src/SimplexCellLists.jl | 10 +++++----- src/collide-forces.jl | 29 +++++++++++++++++------------ test/test-collide-forces.jl | 9 +++++++++ 4 files changed, 31 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index 1880fcf..76e8b5f 100644 --- a/Project.toml +++ b/Project.toml @@ -8,13 +8,11 @@ projects = ["test", "benchmark"] [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" -ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "2" -ChunkSplitters = "3" LinearAlgebra = "1" StaticArrays = "1" julia = "1.12" diff --git a/src/SimplexCellLists.jl b/src/SimplexCellLists.jl index 9fe4139..5b3d899 100644 --- a/src/SimplexCellLists.jl +++ b/src/SimplexCellLists.jl @@ -3,7 +3,6 @@ module SimplexCellLists using LinearAlgebra: LinearAlgebra, cross, dot, ⋅ using StaticArrays: StaticArrays, @SVector, SA, SVector using ArgCheck: @argcheck -using ChunkSplitters: chunks const Vec3 = SVector{3} const Simplex{N,T} = SVector{N, Vec3{T}} @@ -100,10 +99,10 @@ export CollidePairs export N_COLLIDE_PAIRS export empty_no_collide_pairs export CollideObjectTypes -export CollisionPolicy -export filter_object -export filter_pair -export mix_params +public CollisionPolicy +public filter_object +public filter_pair +public mix_params export DefaultObjectParams export DefaultPairParams export DefaultCollisionPolicy @@ -130,6 +129,7 @@ end include("collide-forces.jl") export collide_forces! +public nl_edge_forces! # Precompile the collide forces for the default policy let NL = NeighborLists{DefaultCollisionPolicy, DefaultPairParams} diff --git a/src/collide-forces.jl b/src/collide-forces.jl index 6d137fd..efea2f3 100644 --- a/src/collide-forces.jl +++ b/src/collide-forces.jl @@ -8,7 +8,7 @@ Methods dispatch on the index part types of `edge` and on `policy`, so a custom """ function nl_edge_forces! end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, PointIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, PointIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T x1 = pos[edge.a.i] x2 = pos[edge.b.i] d = map(calc_type, x1 - x2) @@ -27,7 +27,7 @@ Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos nothing end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T x1 = pos[edge.a.i] y1 = pos[edge.b.i] y2 = pos[edge.b.i + UInt32(1)] @@ -51,7 +51,7 @@ Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos nothing end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T x1 = pos[edge.a.i] y1 = pos[edge.b.i] y2 = pos[edge.b.j] @@ -85,7 +85,7 @@ Returns (do_force::Bool, u, fp1, fq0, fq1) where: - fq1: force on Q1 (without k factor) Force on P0 can be computed as -(fp1+fq0+fq1) """ -@inline function line_line_force_core(P0, P1, Q0, Q1, calc_type, switchover_scale_unitless, L0) +@inline function line_line_force_core(P0, P1, Q0, Q1, calc_type::T, switchover_scale_unitless, L0) where T P = map(calc_type, P1-P0) Q = map(calc_type, Q1-Q0) P0mQ0 = map(calc_type, P0-Q0) @@ -196,7 +196,7 @@ Force on P0 can be computed as -(fp1+fq0+fq1) end end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, CLineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T P0 = pos[edge.a.i] P1 = pos[edge.a.i + UInt32(1)] Q0 = pos[edge.b.i] @@ -216,7 +216,7 @@ Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos nothing end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{CLineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T P0 = pos[edge.a.i] P1 = pos[edge.a.i + UInt32(1)] Q0 = pos[edge.b.i] @@ -236,7 +236,7 @@ Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos nothing end -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{LineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{LineIdxPart, LineIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T P0 = pos[edge.a.i] P1 = pos[edge.a.j] Q0 = pos[edge.b.i] @@ -265,7 +265,7 @@ Based on https://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pd // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt // Version: 6.0.2022.01.06 =# -Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, TriangleIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type) +Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos, edge::NeighborListEdge{PointIdxPart, TriangleIdxPart, DefaultPairParams}, policy::DefaultCollisionPolicy, calc_type::T) where T P = pos[edge.a.i] B = pos[edge.b.i] E0 = map(calc_type, pos[edge.b.j] - B) @@ -335,10 +335,15 @@ Base.@propagate_inbounds function nl_edge_forces!(force_energy::ForceEnergy, pos end function nl_forces!(force_energy::ForceEnergy, pos, nl, policy::CollisionPolicy, calc_type::T, chunk, nthreads) where T - chunk > length(nl) && return - for edge in chunks(nl; n=nthreads)[chunk] - nl_edge_forces!(force_energy, pos, edge, policy, calc_type) + # Split 1:length(nl) into nthreads contiguous chunks with sizes differing + # by at most one; chunks past the end are empty. + q, r = divrem(length(nl), nthreads) + start = (chunk - 1)*q + min(chunk - 1, r) + 1 + stop = chunk*q + min(chunk, r) + for i in start:stop + nl_edge_forces!(force_energy, pos, nl[i], policy, calc_type) end + nothing end """ @@ -356,7 +361,7 @@ its own `force_energy`, and the per-thread accumulators are merged afterwards with [`combine_force_energy!`](@ref). Every edge is visited exactly once across the chunks. """ -function collide_forces!(force_energy::ForceEnergy, pos, s::NeighborLists, calc_type; chunk=1, nthreads=1) +function collide_forces!(force_energy::ForceEnergy, pos, s::NeighborLists, calc_type::T; chunk=1, nthreads=1) where T nl_forces!(force_energy, pos, s.PPNL, s.policy, calc_type, chunk, nthreads) nl_forces!(force_energy, pos, s.PCNL, s.policy, calc_type, chunk, nthreads) nl_forces!(force_energy, pos, s.PLNL, s.policy, calc_type, chunk, nthreads) diff --git a/test/test-collide-forces.jl b/test/test-collide-forces.jl index 02a6132..d7424c9 100644 --- a/test/test-collide-forces.jl +++ b/test/test-collide-forces.jl @@ -228,5 +228,14 @@ end @test get_energy(fe_chunked) ≈ get_energy(fe) @test get_force(fe_chunked) ≈ get_force(fe) end + + # empty neighbor lists are fine for any chunk + empty_nl = NeighborLists(policy) + fe_empty = ForceEnergyFloat64(length(positions)) + for chunk in 1:3 + collide_forces!(fe_empty, positions, empty_nl, Float64; chunk, nthreads=3) + end + @test iszero(get_energy(fe_empty)) + @test all(iszero, get_force(fe_empty)) end nothing