2D Material Point Method
2D 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:
The following papers are implemented here:
The following libraries are includes in the ext/
directory:
This project additionally requires the following libraries:
Eigen
can be used to replace Algebra
, but it is not as fast. The source code is in ext/Eigen/MPM2D/src/
(not updated).The followings are optional dependencies :
The code, located in src/
, is structured as following:
main.cpp
: OpenGL context. Run simulation.solver.h
and solver.cpp
: MPM algorithm functions (transfers and updates). Rendering and WriteToFile.node.h
and node.cpp
: Class for grid nodes.border.h
and border.cpp
: Class for 2D linear borders. Collision and Friction.particle.h
and particle.cpp
: Class and subclasses for particles and materials. Constitutive model and deformation functions.constants.h
: Option control and global constants.
Here are the main features of this implementation:
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:
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;
}
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
}
The shape of the domain can be changed, but is has to follow this rules:
CUB
; X_GRID - CUB
] x [CUB
; Y_GRID - CUB
], where CUB
is the range of the interpolation function (2 for Cubic, 1.5 for Quadratic).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;
}
Here is a list of different options available. They can be modify in the constants.h
file.
// Size of the domain
const static int X_GRID = 128;
const static int Y_GRID = 64;
// Select Particle subclass (material type). [Water], [DrySand], [Snow], [Elastic]
#define Material NewMaterial
// Interpolation type: [1] Cubic - [2] Quadratic
#define INTERPOLATION 1
// Time-step (typically about 1e-4)
const static float DT = 0.0001f;
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)