Add collide forces

Project.toml

@@ -1,7 +1,10 @@
 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"
@@ -13,6 +16,3 @@ ArgCheck = "2"
 LinearAlgebra = "1"
 StaticArrays = "1"
 julia = "1.12"
-
-[workspace]
-projects = ["test", "benchmark"]

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, 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.
+- 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 +34,50 @@ 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. 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
+
+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.

src/SimplexCellLists.jl

@@ -99,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
@@ -127,4 +127,22 @@ let NL = NeighborLists{DefaultCollisionPolicy, DefaultPairParams},
     end
 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}
+    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