Deformable Particles

This guide covers the DeformableParticleModel, a bonded-force model that turns a collection of point particles (vertices) into elastic, deformable bodies.

A deformable particle is defined by its mesh: a set of vertices connected by elements (triangles in 3D, segments in 2D). The model computes elastic forces from the mesh geometry. It penalizes deviations in element measure (area/length), body content (volume/area), bending angle, edge length, and surface tension.

This guide shows how to create, configure, and extend deformable particles.

Creating a Deformable Particle

Create a deformable particle with create(). Pass the registered name "deformable_particle_model" and the mesh topology.

Vertices define particle positions, and elements define connectivity. If you do not provide reference (stress-free) quantities, the constructor computes them from vertices. If you provide all required reference quantities, the vertices argument is optional.

The connectivity arrays of the deformable particle contain the particles’ array indices from State.

For the simulation to run correctly, all vertices the deformable model references must exist in the simulation state.

import jax
import jax.numpy as jnp
import jaxdem as jdem

# A simple square boundary in 2D: 4 vertices connected by 4 segments.
vertices_2d = jnp.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]], dtype=float)
elements_2d = jnp.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=int)
edges_2d = elements_2d  # In 2D the edges often coincide with elements.
adjacency_2d = jnp.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=int)

dp = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    edges=edges_2d,
    element_adjacency=adjacency_2d,
    em=1.0,
    eb=0.5,
    el=0.3,
    gamma=0.1,
)
print("Created DP:", type(dp).__name__)

Created DP: DeformableParticleModel

Passing the Model to the System

There are two equivalent ways to attach a deformable particle model to a System.

When you provide a bonded-force model, create() automatically registers its force and energy functions with the ForceManager. The force manager then computes the bonded forces alongside any other forces (contact, gravity, custom) during each time step. You do not pass them manually to the force manager.

Option 1 — pass the model object directly:

state = jdem.State.create(pos=vertices_2d)
system = jdem.System.create(state.shape, bonded_force_model=dp)

Option 2 — pass the registered type name and keyword arguments. create() will build the model internally:

system = jdem.System.create(
    state.shape,
    bonded_force_model_type="deformable_particle_model",
    bonded_force_model_kw={
        "vertices": vertices_2d,
        "elements": elements_2d,
        "edges": edges_2d,
        "element_adjacency": adjacency_2d,
        "em": 1.0,
        "eb": 0.5,
        "el": 0.3,
        "gamma": 0.1,
    },
)

If you pass both, the model object takes precedence and create() ignores the keyword arguments.

Coefficient Broadcasting

Pass every coefficient as a scalar or as a full array. The constructor broadcasts scalar values to the shape of the corresponding geometric entity:

Coefficient	Target shape	Description
em	(M,) — per element	Measure (area/length) stiffness
gamma	(M,) — per element	Surface/line tension
eb	(A,) — per adjacency pair	Bending stiffness
el	(E,) — per edge	Edge length stiffness
ec	(K,) — per body	Content (volume/area) stiffness

ec is special: it is a per-body coefficient, not per-element. The elements_id array maps each element to its parent body, so the model knows which ec value to read for each element’s content contribution. When only one body is present and you do not provide elements_id, ec must have shape (1,).

dp_scalar = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    edges=edges_2d,
    em=2.0,  # broadcast to shape (4,)
    el=0.3,  # broadcast to shape (4,)
    gamma=0.1,  # broadcast to shape (4,)
)
print("em shape:", dp_scalar.em.shape)  # (4,)
print("el shape:", dp_scalar.el.shape)  # (4,)
print("gamma shape:", dp_scalar.gamma.shape)  # (4,)

em shape: (4,)
el shape: (4,)
gamma shape: (4,)

Lazy Array Creation

The constructor only allocates the arrays it needs. If a coefficient is None (i.e. not provided), the constructor does not store the corresponding topology and reference arrays, even if you passed them. This keeps the model small when only some energy terms are active.

dp_edges_only = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    edges=edges_2d,
    element_adjacency=adjacency_2d,
    el=0.3,  # Only edge springs are active.
)
print("elements stored?", dp_edges_only.elements is not None)  # False
print("edges stored?", dp_edges_only.edges is not None)  # True
print("adjacency stored?", dp_edges_only.element_adjacency is not None)  # False

elements stored? False
edges stored? True
adjacency stored? False

Here the constructor discards elements and element_adjacency, because no coefficient that needs them (em, ec, gamma, eb) was provided.

2D vs 3D Differences

JaxDEM supports both 2D and 3D deformable particles. The key differences are:

Concept	2D	3D
Elements	Segments (M, 2)	Triangles (M, 3)
Measure	Segment length	Triangle area
Content	Enclosed area	Enclosed volume
Bending	Angle at shared vertex	Dihedral angle at shared edge
element_adjacency_edges	Not needed (automatically inferred)	Optional (A, 2) — vertex IDs of the shared edge (auto-inferred if not provided)

The constructor infers the dimension from the vertices. vertices.shape[-1] sets whether the model operates in 2D or 3D, and this must be consistent with elements.shape[-1].

In 3D, each adjacency pair shares an edge (two vertices), and element_adjacency_edges stores those vertex IDs. If you do not provide them, the constructor infers them from the element connectivity. In 2D, adjacencies share a single vertex and the edges array is not needed.

# A minimal 3D example: a tetrahedron with 4 triangular faces.
vertices_3d = jnp.array(
    [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 1.0, 0.0], [0.5, 0.5, 1.0]],
    dtype=float,
)
elements_3d = jnp.array(
    [[0, 1, 2], [0, 1, 3], [1, 2, 3], [0, 2, 3]],
    dtype=int,
)
adjacency_3d = jnp.array(
    [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]],
    dtype=int,
)

dp_3d = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_3d,
    elements=elements_3d,
    element_adjacency=adjacency_3d,
    em=1.0,
    eb=0.5,
)
print("3D elements shape:", dp_3d.elements.shape)  # (4, 3)
print("3D adjacency_edges shape:", dp_3d.element_adjacency_edges.shape)  # (6, 2)

3D elements shape: (4, 3)
3D adjacency_edges shape: (6, 2)

The elements_id Field

When a single deformable particle model contains multiple bodies, the elements_id array identifies which body each element belongs to. The content energy term (ec) requires this mapping, because it is a per-body quantity: the model sums the partial content contribution of each element per body using elements_id.

elements_id has shape (M,) and contains integer body indices. For example, if you have two bodies with 3 and 2 elements respectively:

elements_id = jnp.array([0, 0, 0, 1, 1])
ec = jnp.array([0.5, 0.8])  # one value per body

Without elements_id, the model treats all elements as one body (body 0), and ec must have shape (1,).

dp_two_bodies = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    elements_id=jnp.array([0, 0, 1, 1]),
    ec=jnp.array([0.5, 0.8]),
)
print("ec shape:", dp_two_bodies.ec.shape)  # (2,)
print("elements_id:", dp_two_bodies.elements_id)

ec shape: (2,)
elements_id: [0 0 1 1]

Connectivity Masking

Contact forces should skip sphere pairs (or deformable-particle vertices) connected by a bond, because the bonded force model already handles their interaction.

The bond_id field in the State filters contacts. For each sphere, bond_id stores the unique IDs of the spheres it connects to. Connected sphere pairs do not collide.

You can specify connections for each sphere as lists of neighbor unique IDs. The constructor automatically symmetrizes these connections and pads the resulting array with -1.

For example, if we have 3 spheres in a line where 0-1 and 1-2 are connected:

state = jdem.State.create(
    pos=jnp.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]]),
    rad=jnp.array([1.1, 1.1, 1.1]),
    bond_id=[[1], [0, 2], [1]],  # uneven lists of connected unique IDs
)
print("Constructed bond_id:\n", state.bond_id)

Constructed bond_id:
 [[ 1 -1]
 [ 0  2]
 [ 1 -1]]

If a sphere has no connections, you can pass None (or empty lists), which pads everything with -1:

state_no_bonds = jdem.State.create(
    pos=jnp.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]]),
    rad=jnp.array([1.1, 1.1, 1.1]),
    bond_id=None,
)
print("No bonds bond_id:\n", state_no_bonds.bond_id)

No bonds bond_id:
 [[-1]
 [-1]
 [-1]]

Adding and Merging Deformable Particles

Just like State, deformable particle models support add and merge operations for building up complex configurations from smaller pieces.

add() creates a new body from raw arrays and merges it into an existing model. It is equivalent to calling create followed by merge.

dp_base = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    edges=edges_2d,
    em=1.0,
    el=0.3,
)

# Add a second body.
new_verts = vertices_2d + jnp.array([2.0, 0.0])
dp_extended = dp_base.add(
    dp_base,
    vertices=new_verts,
    elements=elements_2d,
    edges=edges_2d,
    em=2.0,
    el=0.5,
)
print("Elements after add:", dp_extended.elements.shape)  # (8, 2)
print("Edges after add:", dp_extended.edges.shape)  # (8, 2)

Elements after add: (8, 2)
Edges after add: (8, 2)

merge() concatenates two existing models. It shifts vertex indices and body IDs so that references stay consistent, and shifts element_id to keep it unique across the merged model. Each merged container therefore represents new bodies with new element IDs. When one side has a term and the other does not, merge pads missing coefficients with 0 and missing reference values with 1.

dp_a = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    em=1.0,
    gamma=0.1,
)

dp_b = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    elements=elements_2d,
    edges=edges_2d,
    el=0.5,
)

dp_merged = dp_a.merge(dp_a, dp_b)
print("Merged em:", dp_merged.em)  # em padded with 0 for dp_b's elements
print("Merged el:", dp_merged.el)  # el padded with 0 for dp_a's edges

Merged em: [1. 1. 1. 1.]
Merged el: [0.5 0.5 0.5 0.5]

Batched Simulations with vmap

Deformable particles work with jax.vmap(), so you can run many independent simulations in parallel. Each simulation gets its own State and System (including its own bonded model).

def create_sim(_i: jax.Array) -> tuple[jdem.State, jdem.System]:
    state = jdem.State.create(pos=vertices_2d)
    dp_model = jdem.BondedForceModel.create(
        "deformable_particle_model",
        vertices=state.pos,
        elements=elements_2d,
        edges=edges_2d,
        element_adjacency=adjacency_2d,
        em=[1.0],
        eb=[0.5],
        el=0.3,
        gamma=0.1,
    )
    system = jdem.System.create(state.shape, bonded_force_model=dp_model)
    return state, system


# Build a batch of 8 independent simulations.
states, systems = jax.vmap(create_sim)(jnp.arange(8))
print("Batched pos shape (B, N, dim):", states.pos.shape)

# Advance all simulations by 5 steps in parallel.
states, systems = systems.step(states, systems, n=5)
print("After stepping:", states.pos.shape)

Batched pos shape (B, N, dim): (8, 4, 2)
After stepping: (8, 4, 2)

A Note on Edges

The edges array and the el coefficient define spring connections between vertex pairs. These are fully independent from the elements connectivity: you can define edge springs that do not correspond to any mesh element. This lets you model springs outside the mesh geometry, such as cross-bracing springs or tethers between non-adjacent vertices.

Because edges are independent, use them alone (without elements) or in combination with any other energy term.

dp_extra_springs = jdem.BondedForceModel.create(
    "deformable_particle_model",
    vertices=vertices_2d,
    edges=jnp.array([[0, 2], [1, 3]]),  # Diagonal springs (not mesh edges).
    el=0.5,
)
print("Diagonal springs — edges:", dp_extra_springs.edges)
print("Diagonal springs — elements stored?", dp_extra_springs.elements is not None)

Diagonal springs — edges: [[0 2]
 [1 3]]
Diagonal springs — elements stored? False

VTK Output

The VTKWriter detects when a DeformableParticleModel is attached to the system and writes additional VTK files:

• deformable_elements — the mesh elements (triangles/segments) with per-cell data: elements_id, ec, gamma, initial_element_measures, current_element_measures, partial_content, and element_normals.

• deformable_edges — the edge springs with per-cell data: initial_edge_lengths, current_edge_lengths, and el.

• deformable_edge_adjacencies — the adjacency pairs (hinge edges in 3D, hinge vertices in 2D) with per-cell data: initial_bendings, current_bendings, and eb.

The writer only writes the files whose energy terms are active. For example, if el is None, it skips the edges writer. You can load this output in ParaView for visualization and debugging.