jaxdem.utils.surface_properties#

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

The target (“central”) clump stays fixed. The sweep places a tracer clump at a sequence of approach directions on a sphere (3D) or circle (2D) surrounding the target. At every approach direction, the sweep pushes the tracer toward the target along the center-to-center direction until the two clumps reach a user-specified geometric overlap – the maximum pairwise sphere overlap delta = r_i + r_j - |x_i - x_j| over all (central-sphere, tracer-sphere) pairs. At the target overlap, the sweep decomposes the interaction force into normal/tangential components with respect to the center-to-center axis, which gives an effective friction coefficient mu = |F_t| / |F_n|.

The sweep repeats 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, sweeping over (facing direction on S^2, roll angle about that facing axis) gives full SO(3) orientation coverage, 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. The function ignores their initial orientations (it resets both 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. The function ignores their initial orientations (it resets both 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, the function builds a default static measurement system (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" the function spaces the points equally on S^1 (2D) and places them on 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). The function pairs each facing 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 get 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 values improve GPU usage. Smaller values cut 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). F is the total force on the central clump, so this is the clump-level coefficient (the same quantity as compute_clump_pair_friction()). Reported as NaN for any probe that fails to establish contact (no normal force). This is deliberately distinct from a genuine frictionless mu == 0.

  • mu_asperity – per-asperity friction coefficient, shape (n_points, *orientation_shape, n_central_spheres). Entry [..., s] decomposes the force on central vertex sphere s alone along the same center-to-center axis used for mu. The last axis follows the sphere order of central_state. NaN where sphere s bears no force (it is not one of the contacting asperities), so np.isnan selects the inactive asperities and np.nanmax(..., axis=-1) gives the worst active asperity per probe. inf in the rare case where the asperity force is exactly tangential to the COM axis, which is a real divergence of the ratio rather than a failed probe.

  • 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, equal to 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"). It records exactly how the function drew the grids, 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.

mu and mu_asperity differ only in the order of summation and decomposition: mu decomposes the summed force, mu_asperity decomposes each asperity force separately. Because the normal axis is shared by every asperity of a probe, the decomposition is linear, and for probes with a single force-bearing asperity (n_central_contacts == 1) the two are identical. With several asperities the tangential components partially cancel in the sum, so mu is bounded above by the normal-force-weighted mean of the active mu_asperity values whenever those asperities all push outward along the COM axis. When some push inward the normal components cancel too and mu can exceed every mu_asperity; those probes are worth inspecting separately.

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).