Simulation Setup

In this section, we define the physical and numerical setup required for implementing a minimal MPM simulation of Two Colliding Elastic Blocks in 2D. We walk through the definition of simulation properties, initialization of particle positions and velocities, and data structures used throughout the simulation.

Physical and Numerical Parameters

We begin by setting up the discretization of the simulation domain and the material parameters of the block:

Implementation 30.1.1 (Physical and Numerical Parameters, simulator.py).

# simulation setup
grid_size = 128 # background Eulerian grid's resolution, in 2D is [128, 128]
dx = 1.0 / grid_size # the domain size is [1m, 1m] in 2D, so dx for each cell is (1/128)m
dt = 2e-4 # time step size in second
ppc = 8 # average particles per cell

density = 1000 # mass density, unit: kg / m^3
E, nu = 1e4, 0.3 # block's Young's modulus and Poisson's ratio
mu, lam = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters

These parameters define a uniform dense background grid, particle resolution, and time integration step size. The entire simulation domain spans from [0, 0] to [1, 1] meters, and we aim for around 8 particles per grid cell on average. The blocks are set to have mass density at $1000 k g / m^{3}$ , Young's modulus at $1 0^{4} P a$ , and Poisson's ratio at $0.3$ .

Initial Particle Sampling and Scene Setup

We sample particles from two rectangular regions using uniform grid sampling. These two boxes are placed symmetrically on the left and right sides of the domain and are initialized with opposite velocities to simulate a head-on collision.

Compared to Poisson disk sampling, uniform sampling is easier to implement for analytic shapes, such as boxes and spheres, due to its structured nature and simple parametrization. However, this regularity can lead to aliasing artifacts, such as visible patterns or striping in the simulation, which may introduce unnatural structured noise into the result.

Here we adopt uniform sampling for simplicity and clarity, keeping the focus on the MPM pipeline itself.

Implementation 30.1.2 (Initial Particle Sampling and Scene Setup, simulator.py).

# uniformly sampling material particles
def uniform_grid(x0, y0, x1, y1, dx):
    xx, yy = np.meshgrid(np.arange(x0, x1 + dx, dx), np.arange(y0, y1 + dx, dx))
    return np.column_stack((xx.ravel(), yy.ravel()))

box1_samples = uniform_grid(0.2, 0.4, 0.4, 0.6, dx / np.sqrt(ppc))
box1_velocities = np.tile(np.array([10.0, 0]), (len(box1_samples), 1))
box2_samples = uniform_grid(0.6, 0.4, 0.8, 0.6, dx / np.sqrt(ppc))
box2_velocities = np.tile(np.array([-10.0, 0]), (len(box1_samples), 1))
all_samples = np.concatenate([box1_samples, box2_samples], axis=0)
all_velocities = np.concatenate([box1_velocities, box2_velocities], axis=0)

Each block consists of uniformly distributed material points representing a homogeneous elastic body. The left block is given an initial velocity of ([+10, 0]) m/s, and the right block ([-10, 0]) m/s, setting up a symmetric, head-on collision scenario with zero net linear momentum. This configuration mimics a controlled impact experiment.

Particle and Grid Data Fields

We define data fields to represent the state of each material point (particle) and background grid node. For particles, this includes position, velocity, volume, mass, and deformation gradient, following Material Particles. For the grid, we define nodal mass and velocity fields using dense arrays, which are sufficient for small-scale simulations. These can be further optimized using sparse grid structures—a direction we leave as future work for interested readers.

Implementation 30.1.3 (Particle and Grid Data Fields, simulator.py).

# material particles data
N_particles = len(all_samples)
x = ti.Vector.field(2, float, N_particles) # the position of particles
x.from_numpy(all_samples)
v = ti.Vector.field(2, float, N_particles) # the velocity of particles
v.from_numpy(all_velocities)
vol = ti.field(float, N_particles)         # the volume of particle
vol.fill(0.2 * 0.4 / N_particles) # get the volume of each particle as V_rest / N_particles
m = ti.field(float, N_particles)           # the mass of particle
m.fill(vol[0] * density)
F = ti.Matrix.field(2, 2, float, N_particles)  # the deformation gradient of particles
F.from_numpy(np.tile(np.eye(2), (N_particles, 1, 1)))

# grid data
grid_m = ti.field(float, (grid_size, grid_size))
grid_v = ti.Vector.field(2, float, (grid_size, grid_size))