MPM2D
2D implementation of the Material Point Method.
Sections
Overview
C++ implementation of the Material Point Method.
The Material Point Method is a numerical technique used to simulate the behavior of continuum materials.
The continuum body is described by a number a Lagrangian elements : the material points.
Kinematic equations are solved on the material points
The material points are surrounded by a background Eulerian grid where the dynamic equations are solved.
It can be summarize in 4 main steps:
- Transfer data from particles de grid nodes
- Update node state (apply forces)
- Transfer data from grid nodes to particles
- Update particles state
Papers implemented
The following papers are implemented here:
- Multi-species simulation of porous sand and water mixtures (partially)
- Drucker-Prager Elastoplasticity for Sand Animation
- An angular momentum conserving affine-particle-in-cell method
- A material point method for snow simulation
Dependencies
The following libraries are includes in the ext/
directory:
- Poisson Generator - Initialize particle position.
- Algebra - Little 2D linear algebra library.
This project additionally requires the following libraries:
- OpenGL and GLFW - Visuals.
- Eigen
-
Eigen
can be used to replaceAlgebra
, but it is not as fast. The source code is inext/Eigen/MPM2D/src/
(not updated).
The followings are optional dependencies :
Code structure
The code, located in src/
, is structured as following:
main.cpp
: OpenGL context. Run simulation.solver.h
andsolver.cpp
: MPM algorithm functions (transfers and updates). Rendering and WriteToFile.node.h
andnode.cpp
: Class for grid nodes.border.h
andborder.cpp
: Class for 2D linear borders. Collision and Friction.particle.h
andparticle.cpp
: Class and subclasses for particles and materials. Constitutive model and deformation functions.constants.h
: Option control and global constants.
Implementation
Characteristics:
Here are the main features of this implementation:
- Sand, Water, Snow and purely elastic simulations already implemented
- 2D.
- Affine-Particle-in-Cell (APIC) transfer type.
- B-Spline Quadratic or Cubic interpolation functions (Quadratic is faster, but not as precise).
- Node forces are updated with an explicit method.
- The domain has to be a convex geometry (for collision detection).
Add material type:
It is easy to add a new type of material. In particle.h
and particle.cpp
, create a new subclasse of Particle
. Beside constructors, the subclass must contain the following functions:
- In
particle.h
:
static std::vector<NewMaterial> InitializeParticles() {
// Define initial particle mass, volume, position, velocity and acceleration
std::vector<NewMaterial> outParticles;
// ...
return outParticles;
}
static std::vector<NewMaterial> AddParticles() {
// Define mass, volume, position, velocity and acceleration of particles to add during the simulation
std::vector<NewMaterial> outParticles;
// ...
return outParticles;
}
- In
particle.cpp
:
void NewMaterial::ConstitutiveModel() {
// Update Ap (pre-update deformation gradient)
}
void NewMaterial::UpdateDeformation(const Matrix2f& T) {
// Update deformation gradient.
// T is the sum of the close node velocity gradients.
// Elasticity, Plasticity functions (return-mapping, hardening) ...
}
void NewMaterial::DrawParticle() {
// OpenGL output of particle points
}
Change domain geometry:
The shape of the domain can be changed, but is has to follow this rules:
- It has to be convex.
- It has to be included in [
CUB
;X_GRID - CUB
] x [CUB
;Y_GRID - CUB
], whereCUB
is the range of the interpolation function (2 for Cubic, 1.5 for Quadratic). - Borders have to be straight lines.
To modify the domain, in border.h
, use the InitializeBorders
static function:
static std::vector<Border> InitializeBorders() {
std::vector<Border> outBorders;
std::vector<Vector2f> Corners;
// New border line
Corners.push_back(Vector2f(X1, Y1)); // First point
Corners.push_back(Vector2f(X2, Y2)); // Second point
// type can be [1](sticky), [2](Separating) or [3](Sliding)
// normal has to be oriented inside the domain and normalized
outBorders.push_back(Border(type, normal, Corners));
Corners.clear();
// Add other border
return outBorders;
}
Options
Here is a list of different options available. They can be modify in the constants.h
file.
- Grid:
// Size of the domain
const static int X_GRID = 128;
const static int Y_GRID = 64;
- Particle:
// Select Particle subclass (material type). [Water], [DrySand], [Snow], [Elastic]
#define Material NewMaterial
- Transfer particles <-> grid:
// Interpolation type: [1] Cubic - [2] Quadratic
#define INTERPOLATION 1
// Time-step (typically about 1e-4)
const static float DT = 0.0001f;
- Output (outputs will be generated in the
out/
directory):
// Generate a .mp4 of the OpenGL window
#define RECORD_VIDEO true
// Generate a .ply file with node coordinates
#define WRITE_TO_FILE false
// Draw nodes (active nodes have a different color)
#define DRAW_NODES false // not recommended (slow)