jaxdem.utils.particle_creation#

Utility functions for creating states of GA particles (rigid/DP).

Functions

build_ga_system(particle_radii, ...[, ...])

Catch-all builder: polydisperse GA/DP particles at a target packing fraction.

build_sphere_system(particle_radii, phi, dim, *)

Catch-all builder for a polydisperse sphere packing at a target phi.

create_dp_container(state, *[, em, ec, eb, ...])

Build a DeformableParticleModel (or a plastic variant) from a DP state.

create_ga_state(N, nv, dim, particle_radius, ...)

Build a State of N geometric-asperity bodies.

create_sphere_state(radii, dim, *[, pos, mass])

Build a State of N simple spheres.

distribute_bodies(state, phi, *[, ...])

Randomly distribute each body's centroid in a box sized for phi.

ga_surface_mask(state, *[, group_by])

Per-body surface mask: True iff the node lies on the convex hull of its body.

jaxdem.utils.particle_creation.create_sphere_state(radii: float | Sequence[float] | ndarray | Array, dim: int, *, pos: Sequence[Sequence[float]] | ndarray | Array | None = None, mass: float | Sequence[float] | ndarray | Array = 1.0) → State[source]#

Build a State of N simple spheres.

No Thomson mesh, no Monte-Carlo union-property integration — this is the lightweight sphere primitive. Each sphere is a single-particle body with its own clump_id. Inertia defaults to the solid-disk (2D) / solid-sphere (3D) formula from the given mass and radius (handled inside State.create).

Parameters:
  • radii – Scalar or length-N per-sphere radii.

  • dim – Spatial dimension (2 or 3).

  • pos – Optional (N, dim) centers. If None, the function stacks all spheres at the origin (handy for a single tracer).

  • mass – Scalar or length-N per-sphere masses. Default 1.0.

jaxdem.utils.particle_creation.create_ga_state(N: int, nv: int, dim: int, particle_radius: float, asperity_radius: float, *, particle_type: str = 'clump', core_type: str = 'hollow', aspect_ratio: float | Sequence[float] | None = None, particle_mass: float = 1.0, n_samples: int = 10000000, seed: int | None = None, mesh_type: str = 'thomson', mesh_kwargs: dict[str, Any] | None = None, asperity_dispersity_type: str = 'constant', asperity_dispersity_kwargs: dict[str, Any] | None = None) → State[source]#

Build a State of N geometric-asperity bodies.

The function places surface asperities on a unit shape (set by mesh_type) and adds an optional central core sphere. The union of the asperities (and, if present, the core) defines the body volume from which the function derives the rigid-body / deformable-body properties. Each member is a filled disk or ball; particle_mass is distributed uniformly over the union, with overlaps counted once. Uncovered gaps or cavities remain empty; their shape follows directly from the supplied members.

Parameters:
  • particle_type ({"clump", "dp"}) – “clump” produces rigid bodies with a computed COM / principal-axis orientation / inertia. “dp” produces deformable particles whose nodes share a common bond_id and each carry an equal share of the body mass and union volume. “sphere”-like bodies are the nv = 1 degenerate case of “clump”.

  • core_type ({"hollow", "solid", "phantom"}, optional) – "hollow" (default) uses the union of the solid surface asperities without adding a central sphere; any uncovered interior remains empty. "solid" includes a central core in both the mass properties and the returned contact geometry. "phantom" includes that core in volume, COM, and inertia, then removes it from the returned contact geometry. A phantom core therefore carries solid-core mass properties; it does not represent an empty interior. Recomputing properties from the returned phantom state would integrate the asperities alone.

  • particle_mass (float) – Prescribed total mass of each body’s material union, including the core when core_type is "solid" or "phantom".

  • n_samples (int) – Numerical integration budget for 3D union properties. Disk unions and analytical 3D cases do not depend on it. Other 3D cases choose boundary or volume quadrature according to overlap density, including simultaneous triple and higher-order overlaps. Check convergence of volume, COM, and world-frame inertia at increasing budgets.

  • mesh_type ({"thomson", "icosphere", "fibonacci", "torus", "helix", "arclength", "faceted"}) –

    How the function generates asperity positions on the unit shape.

    • "thomson" (default): generalized Thomson problem on a hyper-ellipsoid (Riesz-energy minimization). Stochastic, near-uniform result seeded by seed.

    • "icosphere": recursive icosahedron subdivision (3D) or regular polygon (2D). Deterministic. nv must be a valid icosphere count.

    • "fibonacci": golden-angle spiral sphere / evenly-spaced circle. Deterministic. Any nv.

    • "torus": 3D genus-1 torus surface.

    • "helix": 3D chiral helix / 2D Archimedean spiral.

    • "arclength": 2D only — equal arc-length spacing on the ellipse/circle perimeter. Closed-form analogue of the converged Thomson ground state in 2D. Deterministic, any nv.

    • "faceted": regular polygon (2D) or icosahedron (3D) with vertex asperities + face/edge-interior fillers. Produces angular (flat-faced, sharp-cornered) particles rather than smooth spheres.

  • mesh_kwargs (dict, optional) – Mesh-specific keyword arguments forwarded to the underlying generator. See each generator in jaxdem.utils.meshes for accepted keys. Examples: {"steps": 5_000, "alpha": 1.0} (thomson — defaults to steps=10_000), {"tube_ratio": 0.25} (torus), {"n_turns": 3, "helix_radius": 0.3} (helix), {"n_facets": 8} (faceted, 2D only).

  • asperity_dispersity_type ({"constant", "bidisperse", "uniform", "truncated_gaussian"}) – How the function draws asperity radii. The function parameterizes every distribution so its theoretical mean is asperity_radius, which keeps core_radius = particle_radius - asperity_radius deterministic across runs and dispersity settings — only the asperity-to-asperity variation changes. With dispersity, the largest asperities naturally protrude beyond particle_radius and the smallest recede. This is by design. See _generate_disperse_asperity_radii() for the per-distribution sampling formulas.

  • asperity_dispersity_kwargs (dict, optional) – Distribution-specific parameters. Defaults: {} (no effect for "constant"), {"size_ratio": 1.4, "fraction_small": 0.5} for "bidisperse", {"size_ratio": 2.0} for "uniform", {"cv": 0.2, "n_sigma": 2.5} for "truncated_gaussian" (where cv = std / mean and cv * n_sigma must be < 1 to keep radii positive).

jaxdem.utils.particle_creation.distribute_bodies(state: State, phi: float, *, domain_type: str = 'periodic', box_aspect: Sequence[float] | None = None, seed: int | None = None, force_tol: float = 1e-12, randomize_orientation: bool = True, group_by: str = 'auto', collider_type: str = 'naive') → tuple[State, Array][source]#

Randomly distribute each body’s centroid in a box sized for phi.

For each body in state (sphere, rigid clump, or deformable particle), the function builds a bounding sphere from its centroid and the max |node - centroid| + rad. It places those bounding spheres uniformly at random in a periodic (or reflective) box sized to the target packing fraction, FIRE-minimizes so nothing overlaps, and translates every body so its centroid lands at the minimized bounding-sphere center. It can also apply a uniformly random per-body rotation.

phi here is the bounding-sphere packing fraction, not the true body packing fraction (which is lower for non-convex clumps / DPs). The target box volume is sum(bounding_sphere_volume) / phi.

Parameters:
  • state – Input state containing one or more bodies. The function infers body grouping from clump_id / bond_id. Use group_by to force one.

  • phi – Target bounding-sphere packing fraction. Must be in (0, 1).

  • domain_type – "periodic" or "reflect" for the analogue sphere minimization domain.

  • box_aspect – Aspect ratios for the box, shape (dim,). Defaults to isotropic.

  • seed – RNG seed for both the initial random centroid placement and the per-body rotation. Drawn randomly if None.

  • force_tol – Absolute force tolerance for bounding-sphere relaxation.

  • randomize_orientation – If True, apply a uniformly-random per-body rotation after placement.

  • group_by – "auto" (default), "clump", or "bond". See _resolve_body_grouping().

  • collider_type – Collider for the analogue sphere minimization, "naive" or "celllist".

Returns:

The input state with per-body translation and (optionally) per-body random orientation applied, and the periodic box size vector.

Return type:

tuple[State, jax.Array]

jaxdem.utils.particle_creation.ga_surface_mask(state: State, *, group_by: str = 'bond') → Array[source]#

Per-body surface mask: True iff the node lies on the convex hull of its body.

For each body, runs scipy.spatial.ConvexHull on the body’s node positions and marks every hull-boundary node as surface, including collinear edge fillers and coplanar face fillers. For bodies with too few nodes to form a hull (<= dim + 1) or that trigger a Qhull numerical failure, the function marks all their nodes as surface.

Exact for convex (or near-convex) bodies: a node is interior iff it lies strictly inside the convex hull of its body. For Thomson-mesh asperity bodies the surface vs. core separation is unambiguous. Non-convex bodies with concave pockets can misclassify pocket nodes as interior (they sit off the global hull).

Parameters:
  • state – State with one or more bodies.

  • group_by – "bond" (default), "clump", or "auto" — see _resolve_body_grouping().

Returns:

(N,) bool mask, True for surface nodes.

Return type:

jax.Array

jaxdem.utils.particle_creation.create_dp_container(state: State, *, em: float | Array | None = None, ec: float | Array | None = None, eb: float | Array | None = None, el: float | Array | None = None, gamma: float | Array | None = None, tau_s: float | Array | None = None, plasticity_type: Literal['edge', 'perimeter', 'bending', 'none', None] = None, group_by: str = 'bond', is_surface: Array | None = None, interior_edges: str | int = 'fan') → Any[source]#

Build a DeformableParticleModel (or a plastic variant) from a DP state.

For each body (grouped by group_by) the function computes the surface topology from the node positions: polar-angle-ordered polygon segments in 2D, or triangulated convex-hull faces in 3D, including coplanar filler nodes. Nonconvex 3D bodies require explicit surface elements through DeformableParticleModel.Create. It wires interior nodes (those flagged False in is_surface) to surface nodes as extra edges carrying the body’s el coefficient. These nodes never participate in surface elements / element_adjacency / initial_bendings.

Energy coefficients (em, ec, eb, el, gamma) are per-body scalars or per-body arrays of shape (n_bodies,), matching the existing generate_ga_deformable_state convention. If you give plasticity_type, you must also give tau_s. The function interprets it per-body (perimeter), per-edge (edge), or per-adjacency (bending) as appropriate.

Parameters:
  • state – State whose nodes make up one or more DPs. Typically produced by create_ga_state(..., particle_type="dp") and optionally placed via distribute_bodies.

  • em – Per-body elastic coefficients. None disables the term.

  • ec – Per-body elastic coefficients. None disables the term.

  • eb – Per-body elastic coefficients. None disables the term.

  • el – Per-body elastic coefficients. None disables the term.

  • gamma – Per-body elastic coefficients. None disables the term.

  • tau_s – Per-body plastic relaxation time. Required whenever plasticity_type is not None.

  • plasticity_type – One of "edge", "perimeter", "bending", or None / "none" for the elastic container.

  • group_by – Body grouping for topology. Default "bond".

  • is_surface – Optional (N,) bool mask of surface nodes. When None (the default) the function computes the mask automatically with ga_surface_mask (per-body convex hull), which is correct for DPs with no core, a solid core, or a phantom core from create_ga_state and for any other approximately-convex body. Pass an explicit mask to control the boundary/interior classification. A mask does not specify nonconvex 3D connectivity.

  • interior_edges – How to connect interior nodes to surface nodes. "fan" (default) connects every interior node to every surface node in its body. An int K connects each interior node to its K nearest surface nodes.

Return type:

DeformableParticleModel or one of its plastic subclasses.

jaxdem.utils.particle_creation.build_ga_system(particle_radii: Sequence[float] | ndarray, vertex_counts: Sequence[int] | ndarray, asperity_radius: float | Sequence[float] | ndarray, phi: float, dim: int, *, particle_type: str = 'clump', core_type: str = 'hollow', aspect_ratio: Any = None, particle_mass: float | Sequence[float] = 1.0, n_property_samples: int = 10000000, mesh_type: str = 'thomson', mesh_kwargs: dict[str, Any] | None = None, asperity_dispersity_type: str = 'constant', asperity_dispersity_kwargs: dict[str, Any] | None = None, domain_type: str = 'periodic', box_aspect: Sequence[float] | None = None, randomize_orientation: bool = True, initial_phi_bb: float = 0.3, compression_step: float = 0.001, compression_force_tol: float = 1e-12, compression_torque_tol: float | None = None, max_n_min_steps_per_outer: int = 200000, compression_progress: bool = False, fire_dt: float = 0.01, dt: float = 0.001, linear_integrator_type: str = 'verlet', rotation_integrator_type: str | None = 'verletspiral', linear_integrator_kw: dict[str, Any] | None = None, rotation_integrator_kw: dict[str, Any] | None = None, force_model_type: str = 'spring', collider_type: str = 'neighborlist', collider_kw: dict[str, Any] | None = None, minimizer: Any = None, minimizer_kw: dict[str, Any] | None = None, target_fn: Any = None, mat_table: Any = None, e_int: float = 1.0, material_type: str = 'elastic', material_kwargs: dict[str, Any] | None = None, matcher_type: str = 'harmonic', matcher_kwargs: dict[str, Any] | None = None, dp_em: float | None = 1.0, dp_ec: float | None = 1.0, dp_eb: float | None = 1.0, dp_el: float | None = 1.0, dp_gamma: float | None = None, dp_tau_s: float | None = None, dp_plasticity_type: Literal['edge', 'perimeter', 'bending', 'none', None] = None, dp_interior_edges: str | int = 'fan', dp_is_surface: Array | None = None, seed: int | None = None) → tuple[Any, ...][source]#

Catch-all builder: polydisperse GA/DP particles at a target packing fraction.

Given per-body polydispersity (radii, vertex counts, aspect ratios, etc.), the builder creates the state and randomly places each body’s bounding sphere at the loose bounding-sphere packing fraction initial_phi_bb. It then energy-minimizes the analogue sphere system, transfers the new centroids back to the bodies, and quasistatically compresses to the target true-body packing fraction phi with a FIRE system. It returns the result wrapped in a System built with the user-requested integrator, collider, and material.

For particle_type="dp" the builder also creates a DeformableParticleModel (or a plastic variant) and wires it into the returned System as the bonded_force_model. In that case it returns (state, system, container). For clumps/spheres it returns (state, system).

Parameters:
  • particle_radii – Per-body arrays (shape (M,)) or scalars for the last two.

  • vertex_counts – Per-body arrays (shape (M,)) or scalars for the last two.

  • asperity_radius – Per-body arrays (shape (M,)) or scalars for the last two.

  • phi – Target true-body packing fraction (uses each body’s union volume, not its bounding sphere).

  • dim – Spatial dimension (2 or 3).

  • particle_type – "clump" (default) or "dp". nv=1 clumps are effectively spheres.

  • core_type – "hollow" (default), "solid", or "phantom". Forwarded to create_ga_state: "hollow" adds no central sphere, while "phantom" retains solid-core mass properties and removes only the core’s contact geometry.

  • aspect_ratio – None (isotropic), scalar, (dim,), (M,), or (M, dim).

  • asperity_dispersity_type – Asperity-radius polydispersity. Forwarded as-is to every per-body create_ga_state call. Theoretical mean of the chosen distribution always equals each body’s asperity_radius, so the bounding/core geometry is unchanged. See create_ga_state for the supported distributions and kwargs.

  • asperity_dispersity_kwargs – Asperity-radius polydispersity. Forwarded as-is to every per-body create_ga_state call. Theoretical mean of the chosen distribution always equals each body’s asperity_radius, so the bounding/core geometry is unchanged. See create_ga_state for the supported distributions and kwargs.

  • domain_type – "periodic" or "reflect" (closed box).

  • initial_phi_bb – Loose bounding-sphere packing fraction used for the initial random placement. It only needs to be below the jamming density of the bounding spheres. The quasistatic compression reaches the target phi afterwards. Defaults to 0.3.

  • compression_step – Packing-fraction increment for quasistatic compression.

  • compression_force_tol – Force and torque tolerances passed to quasistatic_compress_to_packing_fraction.

  • compression_torque_tol – Force and torque tolerances passed to quasistatic_compress_to_packing_fraction.

  • max_n_min_steps_per_outer – Passed to quasistatic_compress_to_packing_fraction.

  • compression_progress – Passed to quasistatic_compress_to_packing_fraction.

  • fire_dt – Time step of the internal FIRE system used for compression. Defaults to 1e-2.

  • dt – Time step for the returned System.

  • linear_integrator_type – Integrator types for the final (returned) System.

  • rotation_integrator_type – Integrator types for the final (returned) System.

  • linear_integrator_kw – Keyword arguments forwarded to the final System’s integrator constructors (e.g. thermostat parameters for "verlet_rescaling").

  • rotation_integrator_kw – Keyword arguments forwarded to the final System’s integrator constructors (e.g. thermostat parameters for "verlet_rescaling").

  • force_model_type – Parameters of the final (returned) System.

  • collider_type – Parameters of the final (returned) System.

  • collider_kw – Parameters of the final (returned) System.

  • mat_table – Optional pre-built material table. If supplied, the function ignores the material and matcher specs.

  • material_type – Material / matcher spec. When mat_table is None the function creates a single material_type with material_kwargs (defaulting to young=e_int, poisson=0.5, density=1.0) and pairs it with the given matcher.

  • material_kwargs – Material / matcher spec. When mat_table is None the function creates a single material_type with material_kwargs (defaulting to young=e_int, poisson=0.5, density=1.0) and pairs it with the given matcher.

  • matcher_type – Material / matcher spec. When mat_table is None the function creates a single material_type with material_kwargs (defaulting to young=e_int, poisson=0.5, density=1.0) and pairs it with the given matcher.

  • matcher_kwargs – Material / matcher spec. When mat_table is None the function creates a single material_type with material_kwargs (defaulting to young=e_int, poisson=0.5, density=1.0) and pairs it with the given matcher.

  • e_int – Material / matcher spec. When mat_table is None the function creates a single material_type with material_kwargs (defaulting to young=e_int, poisson=0.5, density=1.0) and pairs it with the given matcher.

  • dp_em – DP energy parameters. Unused for particle_type="clump".

  • dp_ec – DP energy parameters. Unused for particle_type="clump".

  • dp_eb – DP energy parameters. Unused for particle_type="clump".

  • dp_el – DP energy parameters. Unused for particle_type="clump".

  • dp_gamma – DP plasticity parameters. Unused for particle_type="clump".

  • dp_tau_s – DP plasticity parameters. Unused for particle_type="clump".

  • dp_plasticity_type – DP plasticity parameters. Unused for particle_type="clump".

  • dp_interior_edges – DP energy / topology parameters. Unused for particle_type="clump".

  • dp_is_surface – DP energy / topology parameters. Unused for particle_type="clump".

Returns:

(state, system) for clumps / spheres, or (state, system, container) for DPs.

Return type:

tuple

Raises:

RuntimeError – If compression or mechanical relaxation fails to reach the target.

jaxdem.utils.particle_creation.build_sphere_system(particle_radii: Sequence[float] | ndarray, phi: float, dim: int, *, particle_mass: float | Sequence[float] = 1.0, domain_type: str = 'periodic', box_aspect: Sequence[float] | None = None, initial_phi: float = 0.3, compression_step: float = 0.001, compression_force_tol: float = 1e-12, compression_torque_tol: float | None = None, max_n_min_steps_per_outer: int = 200000, compression_progress: bool = False, fire_dt: float = 0.01, dt: float = 0.001, linear_integrator_type: str = 'verlet', rotation_integrator_type: str = '', linear_integrator_kw: dict[str, Any] | None = None, rotation_integrator_kw: dict[str, Any] | None = None, force_model_type: str = 'spring', collider_type: str = 'neighborlist', collider_kw: dict[str, Any] | None = None, mat_table: Any = None, e_int: float = 1.0, material_type: str = 'elastic', material_kwargs: dict[str, Any] | None = None, matcher_type: str = 'harmonic', matcher_kwargs: dict[str, Any] | None = None, seed: int | None = None, minimizer: Any = None, minimizer_kw: dict[str, Any] | None = None, target_fn: Any = None) → tuple[State, Any][source]#

Catch-all builder for a polydisperse sphere packing at a target phi.

Sphere counterpart to build_ga_system. The builder draws random positions loose (initial_phi — default 0.3 — or the target, whichever is smaller) via random_sphere_configuration, then quasistatically compresses to phi under the same FIRE/spring setup as build_ga_system (the FIRE time step is fire_dt). It returns (state, system) built with the user-requested integrator, collider, and material. The builder forwards linear_integrator_kw / rotation_integrator_kw to the final integrator constructors (e.g. thermostat parameters for linear_integrator_type="verlet_rescaling").

Parameters mirror build_ga_system (minus all the GA/DP-specific knobs that do not apply to bare spheres).

Raises:

RuntimeError – If compression or mechanical relaxation fails to reach the target.