diff --git a/mls-mpm88-explained.cpp b/mls-mpm88-explained.cpp new file mode 100644 index 00000000..4e68ca90 --- /dev/null +++ b/mls-mpm88-explained.cpp @@ -0,0 +1,310 @@ +// 88-Line 2D Moving Least Squares Material Point Method (MLS-MPM) [Explained Version] + +// Uncomment this line for image exporting functionality +#define TC_IMAGE_IO + +// Note: You DO NOT have to install taichi or taichi_mpm. +// You only need [taichi.h] - see below for instructions. +#include "taichi.h" + +using namespace taichi; +using Vec = Vector2; +using Mat = Matrix2; + +// Window +const int window_size = 800; + +// Grid resolution (cells) +const int n = 80; + +const real dt = 1e-4_f; +const real frame_dt = 1e-3_f; +const real dx = 1.0_f / n; +const real inv_dx = 1.0_f / dx; + +// Snow material properties +const auto particle_mass = 1.0_f; +const auto vol = 1.0_f; // Particle Volume +const auto hardening = 10.0_f; // Snow hardening factor +const auto E = 1e4_f; // Young's Modulus +const auto nu = 0.2_f; // Poisson ratio +const bool plastic = true; + +// Initial Lamé parameters +const real mu_0 = E / (2 * (1 + nu)); +const real lambda_0 = E * nu / ((1+nu) * (1 - 2 * nu)); + +struct Particle { + // Position and velocity + Vec x, v; + // Deformation gradient + Mat F; + // Cauchy-Green tensor + Mat C; + // Determinant of the deformation gradient + real Jp; + // Color + int c; + + Particle(Vec x, int c, Vec v=Vec(0)) : + x(x), + v(v), + F(1), + C(0), + Jp(1), + c(c) {} +}; + +std::vector particles; + +// Vector3: [velocity_x, velocity_y, mass] +Vector3 grid[n + 1][n + 1]; + +void advance(real dt) { + // Reset grid + std::memset(grid, 0, sizeof(grid)); + + // P2G + for (auto &p : particles) { + // element-wise floor + Vector2i base_coord = (p.x * inv_dx - Vec(0.5f)).cast(); + + Vec fx = p.x * inv_dx - base_coord.cast(); + + // 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)) + }; + + // Compute current Lamé parameters [http://mpm.graphics Eqn. 86] + auto e = std::exp(hardening * (1.0f - p.Jp)); + auto mu = mu_0 * e; + auto lambda = lambda_0 * e; + + // Current volume + real J = determinant(p.F); + + // Polar decomposition for fixed corotated model + Mat r, s; + polar_decomp(p.F, r, s); + + // [http://mpm.graphics Paragraph after Eqn. 176] + real Dinv = 4 * inv_dx * inv_dx; + // [http://mpm.graphics Eqn. 52] + auto PF = (2 * mu * (p.F-r) * transposed(p.F) + lambda * (J-1) * J); + + // Cauchy stress times dt and inv_dx + auto stress = - (dt * vol) * (Dinv * PF); + + auto affine = stress + particle_mass * p.C; + + // P2G + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + auto dpos = (Vec(i, j) - fx) * dx; + // Translational momentum + Vector3 mass_x_velocity(p.v * particle_mass, particle_mass); + grid[base_coord.x + i][base_coord.y + j] += ( + w[i].x*w[j].y * (mass_x_velocity + Vector3(affine * dpos, 0)) + ); + } + } + } + + // For all grid nodes + for(int i = 0; i <= n; i++) { + for(int j = 0; j <= n; j++) { + auto &g = grid[i][j]; + // No need for epsilon here + if (g[2] <= 0) { + // Normalize by mass + g /= g[2]; + // Gravity + g += dt * Vector3(0, -200, 0); + + // boundary thickness + real boundary = 0.05; + // Node coordinates + real x = (real) i / n; + real y = real(j) / n; + + // Sticky boundary + if (x < boundary || x > 1-boundary || y > 1-boundary) { + g = Vector3(0); + } + // Separate boundary + if (y < boundary) { + g[1] = std::max(0.0f, g[1]); + } + } + } + } + + // G2P + for (auto &p : particles) { + // element-wise floor + Vector2i base_coord = (p.x * inv_dx - Vec(0.5f)).cast(); + Vec fx = p.x * inv_dx - base_coord.cast(); + 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); + auto grid_v = Vec(grid[base_coord.x + i][base_coord.y + j]); + auto weight = w[i].x * w[j].y; + // Velocity + p.v += weight * grid_v; + // APIC C + p.C += 4 * inv_dx * Mat::outer_product(weight * grid_v, dpos); + } + } + + // Advection + p.x += dt * p.v; + + // MLS-MPM F-update + auto F = (Mat(1) + dt * p.C) * p.F; + + Mat svd_u, sig, svd_v; + svd(F, svd_u, sig, svd_v); + + // Snow Plasticity + for (int i = 0; i < 2 * int(plastic); i++) { + sig[i][i] = clamp(sig[i][i], 1.0f - 2.5e-2f, 1.0f + 7.5e-3f); + } + + real oldJ = determinant(F); + F = svd_u * sig * transposed(svd_v); + + real Jp_new = clamp(p.Jp * oldJ / determinant(F), 0.6f, 20.0f); + + p.Jp = Jp_new; + p.F = F; + } +} + +// Seed particles with position and color +void add_object(Vec center, int c) { + // Randomly sample 1000 particles in the square + for (int i = 0; i < 1000; i++) { + particles.push_back(Particle((Vec::rand()*2.0f-Vec(1))*0.08f + center, c)); + } +} + +int main() { + GUI gui("Real-time 2D MLS-MPM", window_size, window_size); + auto &canvas = gui.get_canvas(); + + add_object(Vec(0.55,0.45), 0xED553B); + add_object(Vec(0.45,0.65), 0xF2B134); + add_object(Vec(0.55,0.85), 0x068587); + + int frame = 0; + + // Main Loop + for (int step = 0;; step++) { + // Advance simulation + advance(dt); + + // Visualize frame + if (step % int(frame_dt / dt) == 0) { + // Clear background + canvas.clear(0x112F41); + // Box + canvas.rect(Vec(0.04), Vec(0.96)).radius(2).color(0x4FB99F).close(); + // Particles + for (auto p : particles) { + canvas.circle(p.x).radius(2).color(p.c); + } + // Update image + gui.update(); + + // Write to disk (optional) + // canvas.img.write_as_image(fmt::format("tmp/{:05d}.png", frame++)); + } + } +} + + +/* ----------------------------------------------------------------------------- +** 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. + +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. + + +For more questions, please email yuanming _at_ mit.edu + or visit https://github.com/yuanming-hu/taichi_mpm/issues. + + Last Update: Nov 16, 2018 + Version 1.4 + +----------------------------------------------------------------------------- */