MPM2D

2D implementation of the Material Point Method.

Sections

Overview
Dependencies
Code structure
Implementation
Options

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:

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 replace Algebra, but it is not as fast. The source code is in ext/Eigen/MPM2D/src/ (not updated).

The followings are optional dependencies :

ffmpeg - Output .mp4 videos.
OpenMP - Parallel computing.

Code structure

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.

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], where CUB 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)

Elias-Gu/MPM2D