High-performance moving least squares material point method (MLS-MPM) solver. (ACM Transactions on Graphics, SIGGRAPH 2018)
MIT License
A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid Body Coupling, ACM Transactions on Graphics (SIGGRAPH 2018).
By Yuanming Hu (MIT CSAIL), Yu Fang (Tsinghua University), Ziheng Ge (University of Science and Technology of China), Ziyin Qu (University of Pennsylvania), Yixin Zhu (UCLA), Andre Pradhana (University of Pennsylvania), Chenfanfu Jiang (University of Pennsylvania).
Welcome to join the Discussion Forum.
Update Nov 2021: with the new Taichi programming language, you can run MLS-MPM on GPU with Python 3 after pip install taichi
Supports Linux, OS X and Windows. Tested on Ubuntu 16.04, Ubuntu 18.04, Arch Linux, MinGW, VS2017, OS X 10.11~10.14.
No need to install taichi
or taichi_mpm
- see the end of code for instructions.
//88-Line 2D Moving Least Squares Material Point Method (MLS-MPM)[with comments]
//#define TC_IMAGE_IO // Uncomment this line for image exporting functionality
#include "taichi.h" // Note: You DO NOT have to install taichi or taichi_mpm.
using namespace taichi;// You only need [taichi.h] - see below for instructions.
const int n = 80 /*grid resolution (cells)*/, window_size = 800;
const real dt = 1e-4_f, frame_dt = 1e-3_f, dx = 1.0_f / n, inv_dx = 1.0_f / dx;
auto particle_mass = 1.0_f, vol = 1.0_f;
auto hardening = 10.0_f, E = 1e4_f, nu = 0.2_f;
real mu_0 = E / (2 * (1 + nu)), lambda_0 = E * nu / ((1+nu) * (1 - 2 * nu));
using Vec = Vector2; using Mat = Matrix2; bool plastic = true;
struct Particle { Vec x, v; Mat F, C; real Jp; int c/*color*/;
Particle(Vec x, int c, Vec v=Vec(0)) : x(x), v(v), F(1), C(0), Jp(1), c(c){}};
std::vector<Particle> particles;
Vector3 grid[n + 1][n + 1]; // velocity + mass, node_res = cell_res + 1
void advance(real dt) {
std::memset(grid, 0, sizeof(grid)); // Reset grid
for (auto &p : particles) { // P2G
Vector2i base_coord=(p.x*inv_dx-Vec(0.5_f)).cast<int>();//element-wise floor
Vec fx = p.x * inv_dx - base_coord.cast<real>();
// Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2]
Vec w[3]{Vec(0.5) * sqr(Vec(1.5) - fx), Vec(0.75) - sqr(fx - Vec(1.0)),
Vec(0.5) * sqr(fx - Vec(0.5))};
auto e = std::exp(hardening * (1.0_f - p.Jp)), mu=mu_0*e, lambda=lambda_0*e;
real J = determinant(p.F); // Current volume
Mat r, s; polar_decomp(p.F, r, s); //Polar decomp. for fixed corotated model
auto stress = // Cauchy stress times dt and inv_dx
-4*inv_dx*inv_dx*dt*vol*(2*mu*(p.F-r) * transposed(p.F)+lambda*(J-1)*J);
auto affine = stress+particle_mass*p.C;
for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) { // Scatter to grid
auto dpos = (Vec(i, j) - fx) * dx;
Vector3 mv(p.v * particle_mass, particle_mass); //translational momentum
grid[base_coord.x + i][base_coord.y + j] +=
w[i].x*w[j].y * (mv + Vector3(affine*dpos, 0));
}
}
for(int i = 0; i <= n; i++) for(int j = 0; j <= n; j++) { //For all grid nodes
auto &g = grid[i][j];
if (g[2] > 0) { // No need for epsilon here
g /= g[2]; // Normalize by mass
g += dt * Vector3(0, -200, 0); // Gravity
real boundary=0.05,x=(real)i/n,y=real(j)/n; //boundary thick.,node coord
if (x < boundary||x > 1-boundary||y > 1-boundary) g=Vector3(0); //Sticky
if (y < boundary) g[1] = std::max(0.0_f, g[1]); //"Separate"
}
}
for (auto &p : particles) { // Grid to particle
Vector2i base_coord=(p.x*inv_dx-Vec(0.5_f)).cast<int>();//element-wise floor
Vec fx = p.x * inv_dx - base_coord.cast<real>();
Vec w[3]{Vec(0.5) * sqr(Vec(1.5) - fx), Vec(0.75) - sqr(fx - Vec(1.0)),
Vec(0.5) * sqr(fx - Vec(0.5))};
p.C = Mat(0); p.v = Vec(0);
for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) {
auto dpos = (Vec(i, j) - fx),
grid_v = Vec(grid[base_coord.x + i][base_coord.y + j]);
auto weight = w[i].x * w[j].y;
p.v += weight * grid_v; // Velocity
p.C += 4 * inv_dx * Mat::outer_product(weight * grid_v, dpos); // APIC C
}
p.x += dt * p.v; // Advection
auto F = (Mat(1) + dt * p.C) * p.F; // MLS-MPM F-update
Mat svd_u, sig, svd_v; svd(F, svd_u, sig, svd_v);
for (int i = 0; i < 2 * int(plastic); i++) // Snow Plasticity
sig[i][i] = clamp(sig[i][i], 1.0_f - 2.5e-2_f, 1.0_f + 7.5e-3_f);
real oldJ = determinant(F); F = svd_u * sig * transposed(svd_v);
real Jp_new = clamp(p.Jp * oldJ / determinant(F), 0.6_f, 20.0_f);
p.Jp = Jp_new; p.F = F;
}
}
void add_object(Vec center, int c) { // Seed particles with position and color
for (int i = 0; i < 500; i++) // Randomly sample 1000 particles in the square
particles.push_back(Particle((Vec::rand()*2.0_f-Vec(1))*0.08_f + center, c));
}
int main() {
GUI gui("Real-time 2D MLS-MPM", window_size, window_size);
add_object(Vec(0.55,0.45), 0xED553B); add_object(Vec(0.45,0.65), 0xF2B134);
add_object(Vec(0.55,0.85), 0x068587); auto &canvas = gui.get_canvas();int f=0;
for (int i = 0;; i++) { // Main Loop
advance(dt); // Advance simulation
if (i % int(frame_dt / dt) == 0) { // Visualize frame
canvas.clear(0x112F41); // Clear background
canvas.rect(Vec(0.04), Vec(0.96)).radius(2).color(0x4FB99F).close();// Box
for(auto p:particles)canvas.circle(p.x).radius(2).color(p.c);//Particles
gui.update(); // Update image
// canvas.img.write_as_image(fmt::format("tmp/{:05d}.png", f++));
}
}
} //----------------------------------------------------------------------------
/* -----------------------------------------------------------------------------
** Reference: A Moving Least Squares Material Point Method with Displacement
Discontinuity and Two-Way Rigid Body Coupling (SIGGRAPH 2018)
By Yuanming Hu (who also wrote this 88-line version), Yu Fang, Ziheng Ge,
Ziyin Qu, Yixin Zhu, Andre Pradhana, Chenfanfu Jiang
** Build Instructions:
Step 1: Download and unzip mls-mpm88.zip (Link: http://bit.ly/mls-mpm88)
Now you should have "mls-mpm88.cpp" and "taichi.h".
Step 2: Compile and run
* Linux:
g++ mls-mpm88.cpp -std=c++14 -g -lX11 -lpthread -O3 -o mls-mpm
./mls-mpm
* Windows (MinGW):
g++ mls-mpm88.cpp -std=c++14 -lgdi32 -lpthread -O3 -o mls-mpm
.\mls-mpm.exe
* Windows (Visual Studio 2017+):
- Create an "Empty Project"
- Use taichi.h as the only header, and mls-mpm88.cpp as the only source
- Change configuration to "Release" and "x64"
- Press F5 to compile and run
* OS X:
g++ mls-mpm88.cpp -std=c++14 -framework Cocoa -lpthread -O3 -o mls-mpm
./mls-mpm
** FAQ:
Q1: What does "1e-4_f" mean?
A1: The same as 1e-4f, a float precision real number.
Q2: What is "real"?
A2: real = float in this file.
Q3: What are the hex numbers like 0xED553B?
A3: They are RGB color values.
The color scheme is borrowed from
https://color.adobe.com/Copy-of-Copy-of-Core-color-theme-11449181/
Q4: How can I get higher-quality?
A4: Change n to 320; Change dt to 1e-5; Change E to 2e4;
Change particle per cube from 500 to 8000 (Ln 72).
After the change the whole animation takes ~3 minutes on my computer.
Q5: How to record the animation?
A5: Uncomment Ln 2 and 85 and create a folder named "tmp".
The frames will be saved to "tmp/XXXXX.png".
To get a video, you can use ffmpeg. If you already have taichi installed,
you can simply go to the "tmp" folder and execute
ti video 60
where 60 stands for 60 FPS. A file named "video.mp4" is what you want.
Q6: How was taichi.h generated?
A6: Please check out my #include <taichi> talk:
http://taichi.graphics/wp-content/uploads/2018/11/include_taichi.pdf
and the generation script:
https://github.com/yuanming-hu/taichi/blob/master/misc/amalgamate.py
You can regenerate it using `ti amal`, if you have taichi installed.
Questions go to yuanming _at_ mit.edu
or https://github.com/yuanming-hu/taichi_mpm/issues.
Last Update: March 6, 2019
Version 1.5
----------------------------------------------------------------------------- */
(This is for installing the high-performance 2D/3D solver including MLS-MPM and CPIC. If you want to play with the 88-line MLS-MPM, you don't have to install anything - see here)
Install taichi (legacy branch)
.
Then, in command line
ti install mpm
and it will install this taichi package automatically.
Support coming later.
Every script under the folder scripts/mls-cpic
is executable with python3
.
taichi/outputs/mpm/
;File
node in Houdini to visualize the bgeo
(particles), obj
(3D meshes), poly
(2D polygons) files.(You only need to specify res
in most cases. The default parameters generally work well.)
All parameters:
Vector<dim, int>
) grid resolution. The length of this vector also specifies the dimensionality of the simulation.real
, default: 1e-4
) delta treal
, default: 1.0 / res[0]
)bool
, default: False
): push particles inside level sets out (turn off when you are using sticky level sets)real
, default: 20000.0
) If things do not separate, use this. Typical value: 20000.0.Vector<dim>
, default: (0, -10, 0)
for 3D, (0, -10)
for 2D)real
, default: 0.01
) You can set to 1 / 24
for real frame rate.int
, default: -1
) Number of threads to use. -1
means maximum threads.int
, default: 1000
) Number of frames to simulate.real
, default: 0
) Penetration penalty. Typical values are 1e3
~ 1e4
.bool
, default: True
) Turn on optimization or not. Turning it off if you need to benchmark the less optimized transfers.string
, default: taichi
will use the current file name)bool
, default: False
) Collide rigid body with level set? (Useful for wine & glass.)real
, default: 0
) RPIC damping value should be between 0 and 1 (inclusive).real
, default: 0
) APIC damping value should be between 0 and 1 (inclusive).bool
, default: False
) Log warning when particles get deletedbool
, default: false
) If true
, output particle attributes other than position
.int
, default: 1000
) If bigger than error, sort particle storage in memory every reorder_interval
substeps.int
, default: 1000
) If bigger than error, sort particle storage in memory every reorder_interval
substeps.rigid
, snow
, jelly
, sand
. For non-rigid type, see Particle Attributes
Vector<3, real>
)bool
, default: True
) Is poisson disk sampling or not? Doesn't support type: rigid
, sdf
bool
, default: True
) Is poisson disk periodic or not? Doesn't support 2Dbool
, default: False
) Is poisson disk sampling from source or not (need to defineframe_update
)? Doesn't support 2Drigid
:Vector<3, real>
, default: (0, 0, 0)
) Let the object rotate along with only this axis. Useful for fans or wheels.bool
, must be explicitly specified) Is thin shell or not?real
, default: 0.0
) Coefficient of restitutionreal
, default: 0.0
) Coefficient of frictionreal
, default: 40
for thin shell, 400
for non-thin shell) Volume/area density.Vector<3, real>
, default: (1, 1, 1)
) rescale the object, bigger or smallerVector<dim, real>
, must be explicitly specified)Vector<dim, real>
, default: (0, 0, 0)
)function(real) => Vector<3, real>
)Vector<1 (2D) or 3 (3D), real>
, default: (0, 0, 0)
) Euler anglesVector<1 (2D) or 3 (3D), real>
, default: (0, 0, 0)
)function(real) => Vector<1 (2D) or 3 (3D), real>
) Takes time, returns Euler anglestc.constant_function13(tc.Vector(0, 0, 0))
real
, default: 0
) damping of linear velocity. Typical value: 1
real
, default: 0
) damping of angular velocity. Typical value: 1
jelly
E
: (real
, default: 1e5
) Young's modulusnu
: (real
, default: 0.3
) Poisson's ratiosnow
hardeing
(real
, default: 10
) Hardening coefficientmu_0
(real
, default: 58333.3
) Lame parameterlambda_0
(real
, default: 38888.9
) Lame parametertheta_c
(real
, default: 2.5e-2
) Critical compressiontheta_s
(real
, default: 7.5e-3
) Critical stretchsand
mu_0
(real
, default: 136038
) Lame parameterlambda_0
(real
, default: 204057
) Lame parameterfriction_angle
(real
, default: 30
)cohesion
(real
, default: 0
)beta
(real
, default: 1
)water
k
: (real
, default: 1e5
) Bulk modulusgamma
: (real
, default: 7
)von_mises
youngs_modulus
: (real
, default: 5e3
) Young's modulus (for elasticity)poisson_ratio
: (real
, default: 0.4
) Poisson's ratio (for elasticity, usually no need to change)yield_stress
: (real
, default:1.0
) Radius of yield surface (for plasticity)scripted_motion_3d.py
.rigid_ground_collision.py
.thin_wheels_fans.py
and the wheel is not turning in the right direction, you can try reverse_vertices=True
.real
, in most cases, instead of float
or double
._f
, so that it will have type real
, instead of float
or double
. Example: 1.5_f
(float
or double
depending on build precision) instead of 1.5
(always double
) or 1.5f
(always float
)taichi
(the main lib, master branch) after updating taichi_mpm
.taichi
is up-to-dateCMake
so that all no source files will be detected0.4
means coeff of friction 0.4
-2.4
means coeff of friction 0.4
with slipSyntax:
object1 = mpm.add_particles(...)
object2 = mpm.add_particles(...)
mpm.add_articulation(type='motor', obj0=object1, obj1=object2, axis=(0, 0, 1), power=0.05)
Rotation: enforce two objects to have the same rotation.
type
: rotation
obj0
, obj1
: two objectswater wheel
examplesDistance: enforce two points on two different object to have constant distance
type
: distance
obj0
, obj1
: two objectsoffset0
, offset1
: (Vector<dim, real>
, default: (0, 0, 0)
) offset of two points to the center of mass to each object, in world spacedistance
(real
, default: initial distance between two poitns) target distancepenalty
(real
, default: 1e5
) corrective penaltycrashing_castle
examplesMotor: enforce object to rotate along an axis on another object, and apply torque
type
: motor
obj0
: the wheel
objectobj1
: the body
objectaxis
: (Vector<dim, real>
) the rotation axis in world spacepower
(real
, default: 0
) torque applied per secondmotor.py
Stepper: enforce object to rotate along an axis on another object at a fixed angular velocity
type
: motor
obj0
: the wheel
objectobj1
: the body
objectaxis
: (Vector<dim, real>
) the rotation axis in world spaceangular_velocity
(real
)pd_source = True
in add_particles
initial_velocity
should be a non-zero vectordelta_t=frame_dt
in add_particles
, which enables the frequency of sampling to be consistent with its initial velocityupdate_frequency
.source_sampling.py
, source_sampling_2d.py
Please cite our paper if you use this code for your research:
@article{hu2018mlsmpmcpic,
title={A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid Body Coupling},
author={Hu, Yuanming and Fang, Yu and Ge, Ziheng and Qu, Ziyin and Zhu, Yixin and Pradhana, Andre and Jiang, Chenfanfu},
journal={ACM Transactions on Graphics (TOG)},
volume={37},
number={4},
pages={150},
year={2018},
publisher={ACM}
}