jaxdem.bonded_forces.plastic_deformable_particle#

Implementation of the plastic deformable particle model.

Classes

PlasticDeformableParticleModel(elements, ...)

This model assumes triangular meshes for 3D bodies and linear segment meshes for 2D bodies.
class jaxdem.bonded_forces.plastic_deformable_particle.PlasticDeformableParticleModel(elements: Array | None, edges: Array | None, element_adjacency: Array | None, element_adjacency_edges: Array | None, elements_id: Array | None, initial_body_contents: Array | None, initial_element_measures: Array | None, initial_edge_lengths: Array | None, initial_bendings: Array | None, w_b: Array | None, em: Array | None, ec: Array | None, eb: Array | None, el: Array | None, gamma: Array | None, tau_s: Array | None)#

Bases: BondedForceModel

This model assumes triangular meshes for 3D bodies and linear segment meshes for 2D bodies. Meshes do not need to be closed, but content-based forces will fail if the mesh is not closed. Vertices should be ordered consistently (e.g., CCW) to define positive normals and bending angles.

The container manages the mesh connectivity (elements, edges, etc.) and reference properties (initial measures, contents, lengths, angles) required to compute forces. It supports both 3D (volumetric bodies bounded by triangles) and 2D (planar bodies bounded by line segments).

The general form of the deformable particle potential energy per particle is:

\[\begin{split}E_K &= E_{K,measure} + E_{K,content} + E_{K,bending} + E_{K,edge} + E_{K,surface} \\ E_{K,measure} &= \sum_{m} \frac{em_m}{2{\mathcal{M}_{m,0}}} \left(\mathcal{M}_m - \mathcal{M}_{m,0}\right)^2 \\ E_{K,surface} &= -\sum_{m} \gamma_m \mathcal{M}_m \\ E_{K,content} &= \frac{e_c}{2{\mathcal{C}_{K,0}}} \left(\mathcal{C}_K - \mathcal{C}_{K,0}\right)^2 \\ E_{K,bending} &= \frac{1}{2} \sum_{a} eb_a wb_{a} \left(\theta_a -\theta_{a,0}\right)^2 \\ E_{K,edge} &= \frac{1}{2} \sum_{e} el_e \left(L_e - L_{e,0}\right)^2\end{split}\]

Plasticity:

The plasticity comes from integrating a spring-dashpot equation for the initial length of each edge, evaluated at every force calculation step:

\[L_{e,0}(t+dt) = L_{e,0}(t) + \frac{1}{\tau_{s,e}} (L_e(t) - L_{e,0}(t)) dt\]

where \(L_{e,0}\) is the initial (reference) edge length (initial_edge_lengths), \(L_e\) is the current edge length, \(\tau_{s,e}\) is the relaxation time (tau_s), and \(dt\) is the simulation time step.

Definitions per Dimension:

  • 3D: Measure ($mathcal{M}$) is Face Area; Content ($mathcal{C}$) is Volume; Elements are Triangles.

  • 2D: Measure ($mathcal{M}$) is Segment Length; Content ($mathcal{C}$) is Enclosed Area; Elements are Segments.

Bending normalization is precomputed from the reference configuration:

  • 3D: \(w_b = l_0 / h_0\), where \(l_0\) is the initial shared-edge (hinge) length and \(h_0\) is the initial distance between adjacent face centroids.

  • 2D: \(w_b = 1/h_0\), with \(h_0 = 0.5 (L_{\text{left},0}+L_{\text{right},0})\), where \(L_{\text{left},0}\) and \(L_{\text{right},0}\) are initial lengths of adjacent segments.

Shapes:

  • K: Number of deformable bodies, which can be defined by their vertices and connectivity (elements, edges, etc.)

  • M: Number of boundary elements (\(K \sum_K m_K\))

  • E: Number of unique edges (\(K \sum_K e_K\))

  • A: Number of element adjacencies (\(K \sum_K a_K\))

All mesh properties of each body are concatenated along axis=0. We do not need to distinguish between bodies to compute forces, except for the content-related term, which requires the body ID for each element to sum up the content contributions per body. The body ID for each element is stored in elements_id, which maps each element to its corresponding body.

Coefficient broadcasting:

  • Scalar coefficients are broadcast to the corresponding geometric entities. em and gamma are broadcast to shape (M,), eb to (A,), and el to (E,).

  • ec is special: its shape is (K,), one value per body, where K is the number of unique IDs in elements_id. elements_id[m] maps each element m to the body index used to read ec.

elements: Array | None#

Array of vertex indices forming the boundary elements. Shape: (M, 3) for 3D (Triangles) or (M, 2) for 2D (Segments). Indices refer to the particle unique_id. Vertices correspond to State.pos.

edges: Array | None#

(E, 2). Each row contains the indices of the two vertices forming the edge. Note: In 2D, the set of edges often overlaps with the set of elements (segments).

Type:

Array of vertex indices forming the edges. Shape

element_adjacency: Array | None#

(A, 2). Each row contains the indices of the two elements sharing a connection.

Type:

Array of element adjacency pairs (for bending/dihedral angles). Shape

element_adjacency_edges: Array | None#

(A, 2). This can be independent from edges because we could have extra edge springs that do not correspond to the mesh connectivity.

Type:

Array of vertex IDs forming the shared edge for each adjacency. Shape

elements_id: Array | None#

(M,). elements_id[i] == k means element i belongs to body k.

Type:

Array of body IDs for each boundary element. Shape

initial_body_contents: Array | None#

(K,). Represents Volume in 3D or Area in 2D.

Type:

Array of reference (stress-free) bulk content for each body. Shape

initial_element_measures: Array | None#

(M,). Represents Area in 3D or Length in 2D.

Type:

Array of reference (stress-free) measures for each element. Shape

initial_edge_lengths: Array | None#

(E,).

Type:

Array of reference (stress-free) lengths for each unique edge. Shape

initial_bendings: Array | None#

(A,). Represents Dihedral Angle in 3D or Vertex Angle in 2D.

Type:

Array of reference (stress-free) bending angles for each adjacency. Shape

w_b: Array | None#

(A,). 3D: hinge_length0 / dual_length0 with dual_length0 between face centers. 2D: 1 / dual_length0 with dual_length0 = 0.5 * (L_left0 + L_right0).

Type:

Precomputed bending normalization coefficient per adjacency. Shape

em: Array | None#

(M,). (Controls Area stiffness in 3D; Length stiffness in 2D).

Type:

Measure elasticity coefficient for each element. Shape

ec: Array | None#

(K,). (Controls Volume stiffness in 3D; Area stiffness in 2D).

Type:

Content elasticity coefficient for each body. Shape

eb: Array | None#

(A,).

Type:

Bending elasticity coefficient for each hinge. Shape

el: Array | None#

(E,).

Type:

Edge length elasticity coefficient for each edge. Shape

gamma: Array | None#

(M,).

Type:

Surface/Line tension coefficient for each element. Shape

tau_s: Array | None#

(E,).

Type:

Plastic relaxation time for each edge. Shape

classmethod Create(*, vertices: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, elements: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, edges: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, element_adjacency: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, element_adjacency_edges: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, elements_id: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, initial_body_contents: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, initial_element_measures: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, initial_edge_lengths: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, initial_bendings: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, em: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, ec: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, eb: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, el: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, gamma: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, tau_s: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None, w_b: Array | ndarray | bool | number | bool | int | float | complex | TypedNdArray | None = None) → Self[source]#
static merge(model1: BondedForceModel, model2: BondedForceModel | Sequence[BondedForceModel]) → BondedForceModel[source]#

Merge two or more bonded-force models into one.

Concatenates topology, reference, and coefficient arrays. Vertex indices and body IDs are shifted automatically so that references remain consistent. When one side carries a term that the other does not, missing coefficients are padded with 0 and missing reference values with 1.

Parameters:
Returns:

A new model containing all bodies from both sides.

Return type:

BondedForceModel

static add(model: BondedForceModel, **kwargs: Any) → BondedForceModel[source]#

Create a new body from raw arrays and merge it into an existing model.

This is a convenience wrapper equivalent to calling the concrete Create constructor followed by merge().

Parameters:

  • model (BondedForceModel) – Existing model to extend.

  • **kwargs – Constructor arguments forwarded to the concrete Create method (e.g. vertices, elements, coefficients, …).

Returns:

The extended model.

Return type:

BondedForceModel

static compute_potential_energy_w_aux(pos: jax.Array, state: State, system: System) → tuple[jax.Array, dict[str, jax.Array]][source]#
static compute_potential_energy(pos: jax.Array, state: State, system: System) → jax.Array[source]#
static compute_forces(pos: jax.Array, state: State, system: System) → tuple[jax.Array, jax.Array][source]#
property force_and_energy_fns: tuple[ForceFunction, EnergyFunction, bool][source]#

Build bonded force/energy callables consumed by the force manager.

Returns:

(force_fn, energy_fn, is_com_force) where:

  • force_fn computes bonded force and torque contributions.

  • energy_fn computes bonded potential-energy contributions.

  • is_com_force indicates where force is applied: True for center-of-mass application, False for contact-point application. This has no effect on spheres.

Return type:

Tuple[ForceFunction, EnergyFunction, bool]