jaxdem.utils.surface_properties#

Monte-Carlo-style sampling of the surface of a clump particle with a tracer clump.

The target (“central”) clump is held fixed. A tracer clump is placed at a sequence of approach directions on a sphere (3D) or circle (2D) surrounding the target. At every approach direction, the tracer is pushed toward the target along the center-to-center direction until the two clumps achieve a user-specified geometric overlap – defined as the maximum pairwise sphere overlap delta = r_i + r_j - |x_i - x_j| over all (central-sphere, tracer-sphere) pairs. Once the target overlap is reached, the interaction force is decomposed into normal/tangential components with respect to the center-to-center axis, yielding an effective friction coefficient mu = |F_t| / |F_n|.

This sweep is repeated over a set of tracer orientations, so the map of mu across the target surface represents the tracer-accessible surface area (SASA-like) along with the contact-anisotropy at every sample point.

In 3D, full SO(3) orientation coverage is obtained by sweeping over (facing direction on S^2, roll angle about that facing axis), which correctly handles asymmetric tracers. In 2D the orientation degree of freedom is a single angle.

Functions

compute_surface_properties(central_state, ...)

Sample surface friction / accessibility over the surface of central_state using tracer_state as a probe.
jaxdem.utils.surface_properties.compute_surface_properties(central_state: State, tracer_state: State, target_overlap: float, *, system: System | None = None, n_points: int = 100, n_orientations: int = 1, n_rolls: int = 1, sampling: str = 'lattice', seed: int = 0, separation_tolerance: float = 1e-10, separation_scale: float = 1.1, batch_size: int = 10000) → dict[str, np.ndarray | jax.Array][source]#

Sample surface friction / accessibility over the surface of central_state using tracer_state as a probe.

Parameters:

  • central_state (State) – Single-clump State instances. Their initial orientations are ignored (both are reset per probe: central to identity, tracer to the swept orientation), so the sampling directions live in the body frame of each clump as encoded by its pos_p.

  • tracer_state (State) – Single-clump State instances. Their initial orientations are ignored (both are reset per probe: central to identity, tracer to the swept orientation), so the sampling directions live in the body frame of each clump as encoded by its pos_p.

  • target_overlap (float) – Desired maximum pairwise sphere overlap (r_i + r_j - |x_i - x_j|) between central and tracer at the reported contact configuration. Must be positive; small values correspond to “just barely indented”. The converged state satisfies this exactly up to separation_tolerance.

  • system (System, optional) – Interaction system used to compute forces (for mu). If None, a default static measurement system is built (spring force, elastic material, naive collider, periodic box large enough for the pair).

  • n_points (int) – Exact number of surface sample points (approach directions). With sampling="lattice" these are equispaced on S^1 (2D) / a Fibonacci golden-spiral lattice on S^2 (3D); with sampling="random" they are iid uniform on the circle/sphere.

  • n_orientations (int) – Number of tracer orientations. In 2D: the count of rotation angles (equispaced, or iid U(0, 2 pi) when sampling="random"). In 3D: the count of facing directions on S^2 (Fibonacci lattice, or iid uniform when random); each facing is paired with every roll below, so total orientations are n_orientations * n_rolls.

  • n_rolls (int) – 3D only: number of rolls about the facing axis (equispaced, or iid U(0, 2 pi) when sampling="random"). Must be 1 in 2D (no roll degree of freedom). For asymmetric tracers set this > 1 to obtain full SO(3) coverage.

  • sampling (str) – "lattice" (default) uses the deterministic equispaced / Fibonacci grids described above – current behavior, unchanged. "random" draws every grid (approach directions, facings, rolls / 2D angles) iid uniform, seeded by seed. Random sampling is incommensurate with clump symmetry, so it removes the sampling-vs-particle aliasing that makes the friction PDF collapse to a few delta spikes for symmetric (uniform-asperity) clumps, and it matches theory that assumes a uniform distribution of orientation angles. Reseed (vary seed) for independent draws / error bars.

  • seed (int) – PRNG seed for sampling="random" (ignored for "lattice"). Fixed by default so random runs are reproducible; pass distinct values for independent samples.

  • separation_tolerance (float) – Bisection convergence tolerance on the tracer center-to-center separation. The converged max(overlap) error shrinks linearly with this value. Must be strictly smaller than target_overlap.

  • separation_scale (float) – Safety factor for the upper bound of the bisection bracket.

  • batch_size (int) – Number of probes per vmap call, over the flat n_points * n_orientations * n_rolls probe grid. Larger is better for GPU utilization; smaller cuts peak memory. With the default 10_000 typical sweeps fit in a single kernel launch.

Returns:

A dictionary of stacked ndarrays:

Common

  • mu – friction coefficient |F_t| / |F_n| per probe, shape (n_points, *orientation_shape). Reported as NaN for any probe that fails to establish contact (no normal force); this is deliberately distinct from a genuine frictionless mu == 0.

  • separation – center-to-center distance at target_overlap, same shape as mu. Also NaN for no-contact probes, so np.isnan(separation) (or mu) flags the missed samples.

  • n_central_contacts – int per probe, same shape as mu; number of central-clump vertex spheres with at least one force-bearing external contact at the bisected configuration.

  • n_tracer_contacts – same, for the tracer clump.

  • central_core_contact – bool per probe, same shape as mu; True if any force-bearing contact involves the central clump’s interior core sphere. Always False when the central clump has no core (e.g. a "hollow" clump).

  • tracer_core_contact – same, for the tracer clump’s core.

  • tracer_quaternions – shape (n_points, *orientation_shape, 4); the composed quaternion actually applied to the tracer at the bisected configuration (approach-direction rotation composed with the per-orientation base).

  • central_position – shape (dim,); COM of the central clump (constant across probes; equals system.domain.box_size / 2). Combined with approach_directions and separation this gives the tracer COM as central_position + separation * approach_dir.

  • approach_directions – surface sample directions, shape (n_points, dim).

  • target_overlap – scalar float; echo of the input.

  • dim – int; dimensionality (2 or 3).

  • sampling – str; echo of the sampling mode used.

  • seed – int or None; the PRNG seed for sampling="random" (None for "lattice"), recording exactly how the grids were drawn for reproducibility.

2D only

  • angle_surface – (n_points,); polar angle of each approach direction.

  • tracer_angles – (n_orientations,); the swept tracer rotation angles. Grid shape: (n_orientations,).

3D only

  • theta_surface, phi_surface – (n_points,) each; spherical coordinates of each approach direction.

  • tracer_facings – (n_orientations, 3); body-frame directions sampled on S^2.

  • tracer_facing_theta, tracer_facing_phi – (n_orientations,) each; spherical coords of the facings.

  • tracer_rolls – (n_rolls,); roll angles about each facing axis. Grid shape: (n_orientations, n_rolls).

Return type:

dict

Notes

The orientation_shape is (n_orientations,) in 2D and (n_orientations, n_rolls) in 3D, so mu[i, ...] gives the full orientation map for approach-direction i and any slice along the leading axis gives the surface map for a fixed orientation.

To reproduce the exact contact configuration of probe (i, j[, k]) given the original central_state and tracer_state:

idx = (i, j, k) if dim == 3 else (i, j)
tracer_com = result["central_position"] + (
    result["separation"][idx] * result["approach_directions"][i]
)
tracer_quat = result["tracer_quaternions"][idx]   # (4,)

Apply tracer_com to the tracer’s state.pos_c and tracer_quat to its state.q, leave the central at result["central_position"], and call system.collider.compute_force(state, system) to get the bisected force network.

With sampling="random" the returned approach_directions are not ordered by angle. Consumers that assume ordering – e.g. the 2D SASA perimeter reconstruction, which connects consecutive samples into a polygon – must sort by angle_surface first (3D SASA via ConvexHull is order-independent and needs no change).