Download this testcase.
Tetris¶
This example illustrates how to define a particle problem in which composite rigid bodies—shaped like tetriminoes—fall under gravity into an open box. The pieces are generated randomly and interact through frictional contact, stacking up in a way reminiscent of the classic Tetris game.
Keywords¶
DEM, Rigid bodies, Random generation
Description¶
Each piece (tetrimino) is composed of four circular particles grouped into a single rigid body. Pieces are inserted at random positions and orientations, and fall under gravity into an open container bounded by three fixed walls. The contact solver handles collisions between the falling bodies and the boundaries, leading to complex packing and rotation dynamics.
import sys, os, shutil, random
import numpy as np
from migflow import scontact
random.seed(0)
Output Directory¶
The output directory is prepared for storing the simulation results.
outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)
Piece Generation¶
Function to generate and insert a tetrimino piece composed of four rigidly connected circular particles.
def add_piece(p: scontact.ParticleProblem, r: float, initial_coord: np.ndarray):
"""Create and insert a tetrimino piece into the MigFlow particle problem.
Parameters
----------
p : scontact.ParticleProblem
The particle problem instance to which the piece is added.
r : float
Particle radius.
initial_coord : np.ndarray
Initial coordinates for the lower-left reference position of the piece.
"""
# Each piece is defined by a binary mask over a 2x4 rectangular grid.
models = [
[[1, 0, 0, 0], [1, 1, 1, 0]],
[[0, 0, 0, 1], [0, 1, 1, 1]],
[[1, 1, 1, 1], [0, 0, 0, 0]],
[[0, 1, 1, 0], [0, 1, 1, 0]],
[[0, 0, 1, 1], [0, 1, 1, 0]],
[[0, 1, 0, 0], [1, 1, 1, 0]],
[[0, 1, 1, 0], [0, 0, 1, 1]],
]
piece = random.choice(models)
omega = 10 * np.pi * (-1 + 2 * random.random())
y, x = np.where(piece)
x = r * (initial_coord[0] + 2 * x)
y = r * (initial_coord[1] + 2 * y)
R = np.repeat(r, 4) # particles radii
coord = np.column_stack([x, y]) # particle coordinates
body = p.add_particle_body(coord, R, np.pi * R**2 * rho)
p.body_omega()[body, 0] = omega # initial body rotation
Physical Parameters¶
g = np.array([0, -9.81]) # gravity [m/s²]
rho = 1000 # particle density [kg/m³]
r = 0.05 # particle radius [m]
h = 2 # box height [m]
w = 2 # box width [m]
x0 = np.array([-3, 40]) # initial insertion position
Particle Problem¶
The 2D particle problem is created and the open box boundaries are defined. The container is defined using three rigid walls (left, bottom, right).
p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)
xs0 = np.array([-w / 2, +h / 2])
xs1 = np.array([-w / 2, -h / 2])
xs2 = np.array([+w / 2, -h / 2])
xs3 = np.array([+w / 2, +h / 2])
edges = [[xs0, xs1], [xs1, xs2], [xs2, xs3]]
nodes = [xs0, xs1, xs2, xs3]
bnd_body = p.add_body((0, 0), 0, 0)
for edge in edges:
p.add_segment_to_body(edge[0], edge[1], bnd_body, "bnd")
for n in nodes:
p.add_particle_to_body(n, 0, bnd_body)
Simulation Parameters¶
outf = 20 # number of iterations between output files
dt = 1e-3 # time step [s]
tEnd = 20 # total simulation time [s]
t = 0 # initial time [s]
i = 0 # iteration counter
Computational Loop¶
Pieces fall into the container under gravity. Every 250 iterations, a new piece is generated and inserted.
while t < tEnd:
if i % outf == 0:
print(f"{i =:4d}, {t:.6g}/{tEnd:.6g}")
p.write_mig(outputdir, t) # write results
forces = g * (np.pi * p.r() ** 2 * rho) # external forces (gravity)
p.iterate(dt, forces, tol=1e-4 * r) # contact solver
t += dt
i += 1
# Add a new piece periodically
if i % 250 == 0:
print("add_piece")
add_piece(p, r, x0)
Plot¶
python3 -m migflow.plot.migplot output --actors particles