High-performance moving least squares material point method (MLS-MPM) solver.

High-Performance MLS-MPM Solver with Cutting and Coupling (CPIC) (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.

News

[Introduction & Demo Video] [Paper] [Supplemental Document] [88-Line MLS-MPM]

[SIGGRAPH 2018 Fast Forward] [PDF Slides] [PDF Slides with Notes]

88-Line Version (MIT License) [Download C++ & Javascript versions]


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.

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 ();//element-wise floor 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))}; 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 ();//element-wise floor 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), 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 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 ----------------------------------------------------------------------------- */">
//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
         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
        
         ();
    
         // 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
         
          ();
    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 
           
             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

          

          ----------------------------------------------------------------------------- */
         
        
       

Installing the High-Performance 3D Solver

(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)

Linux and OSX

Install taichi (legacy branch). Then, in command line

ti install mpm

and it will install this taichi package automatically.

Windows

Support coming later.

Run demos

Every script under the folder scripts/mls-cpic is executable with python3.

Visualize the results

  • Outputs are in taichi/outputs/mpm/;
  • Install Houdini Apprentice (which is free);
  • Create a File node in Houdini to visualize the bgeo (particles), obj (3D meshes), poly (2D polygons) files.

Python 3 API

MPM.initialize

(You only need to specify res in most cases. The default parameters generally work well.) All parameters:

  • res: (Vector ) grid resolution. The length of this vector also specifies the dimensionality of the simulation.
  • base_delta_t : (real, default: 1e-4) delta t
  • delta_x: (real, default: 1.0 / res[0])
  • particle_collision (bool, default: False): push particles inside level sets out (turn off when you are using sticky level sets)
  • pushing_force: (real, default: 20000.0) If things do not separate, use this. Typical value: 20000.0.
  • gravity (Vector , default: (0, -10, 0) for 3D, (0, -10) for 2D)
  • frame_dt: (real, default: 0.01) You can set to 1 / 24 for real frame rate.
  • num_threads: (int, default: -1) Number of threads to use. -1 means maximum threads.
  • num_frames: (int, default: 1000) Number of frames to simulate.
  • penalty: (real, default: 0) Penetration penalty. Typical values are 1e3 ~ 1e4.
  • optimized: (bool, default: True) Turn on optimization or not. Turning it off if you need to benchmark the less optimized transfers.
  • task_id: (string, default: taichi will use the current file name)
  • rigid_body_levelset_collision: (bool, default: False) Collide rigid body with level set? (Useful for wine & glass.)
  • rpic_damping: (real, default: 0) RPIC damping value should be between 0 and 1 (inclusive).
  • apic_damping: (real, default: 0) APIC damping value should be between 0 and 1 (inclusive).
  • warn_particle_deletion: (bool, default: False) Log warning when particles get deleted
  • verbose_bgeo: (bool, default: false) If true, output particle attributes other than position.
  • reorder_interval: (int, default: 1000) If bigger than error, sort particle storage in memory every reorder_interval substeps.
  • clean_boundary: (int, default: 1000) If bigger than error, sort particle storage in memory every reorder_interval substeps.
  • ...

MPM.add_particles

  • type: rigid, snow, jelly, sand. For non-rigid type, see Particle Attributes
  • color: (Vector<3, real>)
  • pd : (bool, default: True) Is poisson disk sampling or not? Doesn't support type: rigid, sdf
  • pd_periodic : (bool, default: True) Is poisson disk periodic or not? Doesn't support 2D
  • pd_source : (bool, default: False) Is poisson disk sampling from source or not (need to defineframe_update)? Doesn't support 2D
  • For type rigid:
  • rotation_axis : (Vector<3, real>, default: (0, 0, 0)) Let the object rotate along with only this axis. Useful for fans or wheels.
  • codimensional : (bool, must be explicitly specified) Is thin shell or not?
  • restitution: (real, default: 0.0) Coefficient of restitution
  • friction: (real, default: 0.0) Coefficient of friction
  • density: (real, default: 40 for thin shell, 400 for non-thin shell) Volume/area density.
  • scale: (Vector<3, real>, default: (1, 1, 1)) rescale the object, bigger or smaller
  • initial_position: (Vector , must be explicitly specified)
  • initial_velocity: (Vector , default: (0, 0, 0))
  • scripted_position: (function(real) => Vector<3, real>)
  • initial_rotation: (Vector<1 (2D) or 3 (3D), real>, default: (0, 0, 0)) Euler angles
  • initial_angular_velocity: (Vector<1 (2D) or 3 (3D), real>, default: (0, 0, 0))
  • scripted_rotation: (function(real) => Vector<1 (2D) or 3 (3D), real>) Takes time, returns Euler angles
  • (Translational/Rotational) static objects are also considered as scripted, but with a fixed scripting function i.e. tc.constant_function13(tc.Vector(0, 0, 0))
  • linear_damping: (real, default: 0) damping of linear velocity. Typical value: 1
  • angular_damping: (real, default: 0) damping of angular velocity. Typical value: 1
  • ...

Particle Attributes

  • jelly
    • E: (real, default: 1e5) Young's modulus
    • nu: (real, default: 0.3) Poisson's ratio
  • snow
    • hardeing (real, default: 10) Hardening coefficient
    • mu_0 (real, default: 58333.3) Lame parameter
    • lambda_0 (real, default: 38888.9) Lame parameter
    • theta_c (real, default: 2.5e-2) Critical compression
    • theta_s (real, default: 7.5e-3) Critical stretch
  • sand
    • mu_0 (real, default: 136038) Lame parameter
    • lambda_0 (real, default: 204057) Lame parameter
    • friction_angle (real, default: 30)
    • cohesion (real, default: 0)
    • beta (real, default: 1)
  • water
    • k: (real, default: 1e5) Bulk modulus
    • gamma: (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)
  • ...

Script Examples

  • Scripted motion: scripted_motion_3d.py.
  • Rigid-ground collison: rigid_ground_collision.py.
  • When you're making an rotating wheel example, e.g. thin_wheels_fans.py and the wheel is not turning in the right direction, you can try reverse_vertices=True.
  • ...

Notes

  • Matrices in taichi are column major. E.g. A[3][1] is the element at row 2 and column 4.
  • All indices, unless explicitly specified, are 0-based.
  • Use real, in most cases, instead of float or double.
  • Float point constants should be suffixed with _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)
  • Always pull taichi (the main lib, master branch) after updating taichi_mpm.
  • When a particles moves too close to the boundary (4-8 dx) it will be deleted.
  • Whenever you can any compile/linking problem:
    • Make sure taichi is up-to-date
    • Invoke CMake so that all no source files will be detected
    • Rebuild
  • ...

Friction Coefficient

  • Separate: positive values, 0.4 means coeff of friction 0.4
  • Sticky: -1
  • Slip: -2
  • Slip with friction: -2.4 means coeff of friction 0.4 with slip

Articulation

Syntax:

    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 objects
    • Use case: blabe and wheel in water wheel examples
  • Distance: enforce two points on two different object to have constant distance

    • type: distance
    • obj0, obj1: two objects
    • offset0, offset1: (Vector , default: (0, 0, 0)) offset of two points to the center of mass to each object, in world space
    • distance (real, default: initial distance between two poitns) target distance
    • penalty (real, default: 1e5) corrective penalty
    • Use case: hammer in crashing_castle examples
  • Motor: enforce object to rotate along an axis on another object, and apply torque

    • type: motor
    • obj0: the wheel object
    • obj1: the body object
    • axis: (Vector ) the rotation axis in world space
    • power (real, default: 0) torque applied per second
    • Use case: wheels for cars, and legs for the robot
    • Example: motor.py
  • Stepper: enforce object to rotate along an axis on another object at a fixed angular velocity

    • type: motor
    • obj0: the wheel object
    • obj1: the body object
    • axis: (Vector ) the rotation axis in world space
    • angular_velocity (real)
    • Use case: Fixed-rotation-speed wheels for cars, and legs for the robot

Source Sampling

  • If you want to source particles continuously from a object, please set pd_source = True in add_particles
  • initial_velocity should be a non-zero vector
  • Remember to also set delta_t=frame_dt in add_particles, which enables the frequency of sampling to be consistent with its initial velocity
  • There might be some artifact due to the effect of gravity. You can reduce that artifact by increasing update_frequency.
  • Example: source_sampling.py, source_sampling_2d.py

Mathematical Comparisons with Traditional MPM

Performance

Bibtex

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}
}
Owner
Yuanming Hu
Creator of Taichi; Co-founder and CEO of Taichi Graphics; Ph.D. in Computer Science (MIT CSAIL)
Yuanming Hu
Comments
  • Compilation error

    Compilation error

    Hi When I run "python3 install.py", I get this:

    clang-6/7 is the only supported compiler for Taichi compiler development

    Can you help me?

    Oscariox

  • plasticity

    plasticity

    hello In the case of the snow related plasticity, I encountered the problem of 1 - tetac and 1 + tetas because I wanted to see how it was obtained for the soil of these values (although the empirical values do independ on the hydrostatic stress).

  • Explanation please: (p.x*inv_dx-Vec(0.5_f)).cast<int>();

    Explanation please: (p.x*inv_dx-Vec(0.5_f)).cast();

    Hi,

    just a quick question, what is base_coord after the following line.

    Vector2i base_coord = (p.x*inv_dx-Vec(0.5_f)).cast<int>();
    

    Thanks a lot

  • Command ti install mpm errors

    Command ti install mpm errors

    I was following the instructions to install the High-Performance 3D Solver however the make file for the command ti install mpm failed. I am running on a bash terminal window for a 2018 Macbook Pro. Output: make[1]: *** [projects/mpm/CMakeFiles/taichi_mpm.dir/all] Error 2 make: *** [all] Error 2 Error: Build failed.

    The errors that I get are: error: cannot initialize a variable of type 'char *' with an rvalue of type 'int' error: static_assert failed due to requirement 'bit::is_power_of_two((int)sizeof(GridState<3>))' error: no matching function for call to 'equal' error: no member named 'block_zbits' in 'SPGrid::SPGrid_Mask<6, 6, 2, 12>'

    These errors begin at build percentage 93%. I installed taichi following the instructions in https://taichi.readthedocs.io/en/latest/installation.html and I have all of the correct dependencies and versions of gcc and python.

    Please let me know if I am missing something. Any feedback appreciated.

    screen shot 2019-03-07 at 11 42 11 am
  • Particles leaking from rigid bodies

    Particles leaking from rigid bodies

    Hi. I wanted to make sure that I'm using a correct grid spacing (dx=dimension/resolution) and/or pd's minimum distance. In the following pics (3D and top views) the setting is as default. However, I'm seeing some leakages from the top rigid surface. Please let me know if there is sth that I can figure it out; ty.

    3dview topview

  • error: static assertion failed:

    error: static assertion failed:

    Hi Yuanming, When I installed mpm solver project using command'ti install mpm', several errors appear a lot of times which I think maybe you could help me to solve. Errors are : /taichi/projects/mpm/src/mpm_fwd.h:116:1: error: static assertion failed: GridState<2> size must be POT static_assert(bit::is_power_of_two(sizeof(GridState<2>)), error: ‘using Vector = struct taichi::VectorND<3, float, (taichi::InstSetExt)0> {aka struct taichi::VectorND<3, float, (taichi::InstSetExt)0>}’ has no member named ‘v’ p.pos.v = mm_fmadd_ps(v, delta_t_vec, p.pos.v); ^ Could you please help me have a look at these problems?

  • Installation issues

    Installation issues

    Trying to install the legacy version, but keep getting this when running install.py.

    /home/blazes/taichi/include/taichi/math/linalg.h:372:11: warning: ‘void* memcpy(void*, const void*, size_t)’ writing to an object of type ‘struct taichi::VectorND<4, float, taichi::InstSetExt::AVX2>’ with no trivial copy-assignment; use copy-assignment or copy-initialization instead [-Wclass-memaccess] 372 | memcpy(this, &o, sizeof(*this)); | ~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~ /home/blazes/taichi/include/taichi/math/linalg.h:154:8: note: ‘struct taichi::VectorND<4, float, taichi::InstSetExt::AVX2>’ declared here 154 | struct VectorND : public VectorNDBase<dim__, T, ISE> { | ^~~~~~~~ make[2]: *** [CMakeFiles/taichi_core.dir/build.make:141: CMakeFiles/taichi_core.dir/include/taichi/dynamics/rigid_body.cpp.o] Error 1 make[1]: *** [CMakeFiles/Makefile2:96: CMakeFiles/taichi_core.dir/all] Error 2 make: *** [Makefile:84: all] Error 2 Error: Build failed. Error: installation failed.

  • Installation page not exist

    Installation page not exist

    Clicking on the link provided by "Install taichi (legacy branch)" leads to a page that no longer exists. Where else can I go to install Taichi?

  • Failed on the

    Failed on the "ti install mpm" step

    Hi yuanming,

    I followed the commands (as below) on both my ubuntu16.04 and ubuntu 18.04 (In virtual machine), but I got the same error at the "ti install mpm" step.

    wget https://raw.githubusercontent.com/yuanming-hu/taichi/legacy/install.py python3 install.py ti build ti install mpm

    I didn't run the command "export TC_USE_DOUBLE=1".

    I got an error on the last step("ti install mpm"), which is as follows:

    error********** error: static assertion failed: GridState<2> size must be POT static_assert(bit::is_power_of_two((int)sizeof(GridState<2>)),


    The full error information is as follows:

    full error information*** [ 92%] Building CXX object projects/mpm/CMakeFiles/taichi_mpm.dir/src/mpm_rigid_body.cpp.o In file included from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm.h:24:0, from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/boundary_particle.h:8, from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/boundary_particle.cpp:6: /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm_fwd.h:115:1: error: static assertion failed: GridState<2> size must be POT static_assert(bit::is_power_of_two((int)sizeof(GridState<2>)), ^~~~~~~~~~~~~ In file included from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm.h:24:0, from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm_rigid_body.cpp:8: /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm_fwd.h:115:1: error: static assertion failed: GridState<2> size must be POT static_assert(bit::is_power_of_two((int)sizeof(GridState<2>)), ^~~~~~~~~~~~~ In file included from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm.h:24:0, from /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm.cpp:19: /home/dyfluid/taichi-mpm/installation/taichi/projects/mpm/src/mpm_fwd.h:115:1: error: static assertion failed: GridState<2> size must be POT static_assert(bit::is_power_of_two((int)sizeof(GridState<2>)), ^~~~~~~~~~~~~ projects/mpm/CMakeFiles/taichi_mpm.dir/build.make:86: recipe for target 'projects/mpm/CMakeFiles/taichi_mpm.dir/src/boundary_particle.cpp.o' failed make[2]: *** [projects/mpm/CMakeFiles/taichi_mpm.dir/src/boundary_particle.cpp.o] Error 1 make[2]: *** Waiting for unfinished jobs.... projects/mpm/CMakeFiles/taichi_mpm.dir/build.make:134: recipe for target 'projects/mpm/CMakeFiles/taichi_mpm.dir/src/mpm_rigid_body.cpp.o' failed make[2]: *** [projects/mpm/CMakeFiles/taichi_mpm.dir/src/mpm_rigid_body.cpp.o] Error 1 projects/mpm/CMakeFiles/taichi_mpm.dir/build.make:110: recipe for target 'projects/mpm/CMakeFiles/taichi_mpm.dir/src/mpm.cpp.o' failed make[2]: *** [projects/mpm/CMakeFiles/taichi_mpm.dir/src/mpm.cpp.o] Error 1 CMakeFiles/Makefile2:147: recipe for target 'projects/mpm/CMakeFiles/taichi_mpm.dir/all' failed make[1]: *** [projects/mpm/CMakeFiles/taichi_mpm.dir/all] Error 2 Makefile:129: recipe for target 'all' failed make: *** [all] Error 2 Error: Build failed.


    I have checked other issues and they couldn't solve my problem. Could you please help me have a look at the problem? Is this problem caused by the virtual machine?

    Thank you very much.

  • how to uninstall taichi (legacy branch)

    how to uninstall taichi (legacy branch)

    After experiencing the same issue as #40, I found all my other conda environment are crashed because when I import taichi, seems it will automatically bind to the error taichi rather than the one I pip installed... How can I uninstall the taichi(legacy branch)?

    image

  • Added mls-mpm88 with additional explanation

    Added mls-mpm88 with additional explanation

    Thanks for providing the mls-mpm88.cpp, it is a good MPM in a nutshell example :)

    Prettified the mls-mpm88 example and added additional information on coefficients/steps, hopefully not too many mistakes there!

  • Condition for a grid node in rigid_transfer.cpp for 2D and 3D case.

    Condition for a grid node in rigid_transfer.cpp for 2D and 3D case.

    Hi, I am trying to understand the code, having a problem in understanding rigid_transfer.cpp. Specifically, the condition for determining whether an grid node is valid for a rigid body element.

    If I am understanding correctly, in 3D case described in line 47 to line 51, coord[0] and coord[1] are position of a grid node projected on a triangular mesh and thus it seems to be reasonable to have the condition; 0 <= coord[0] && 0 <= coord[1] && coord[0] + coord[1] <= 1. image

    In 2D case described in line 42 to line 46, however, although the code is applying similar condition to the projected grid point on a elemental bar, the condition seems to be relaxed by 2%. -0.02 <= coord[0] && coord[0] <= 1.02 image

    Could anyone tell me why the difference between 2D and 3D are necessary and where do the numbers -0.02 and 1.02 come from? Thank you so much.

  • Modeling fluid with different viscosity using Taichi_mpm

    Modeling fluid with different viscosity using Taichi_mpm

    Hi,

    I'm a newcomer to computer graphics, so please correct me if I make any mistakes. :)

    I'm trying to model liquid with different viscosity for my new project. I noticed that from Python API end, the two parameters user can change are bulk modulus and gamma, but I found that in this file \src\particle.cpp, line 47 to 49 defined a class ViscoParticle with the following parameters

     real visco_tau;
      real visco_nu;
      real visco_kappa;
    

    Is this what I should be looking for if I want to change the viscosity of fluid? If so, how should I implement this in Taichi Python API? Or will that constrain me to use c++ API only? Many thanks in advance for answering my questions! :))

  • How to open demo results

    How to open demo results

    Hello mpm community, I could successful run banana demo, but have some problems with output. It generated ouputs/mpm/banana folder with .obj, .bgeo files (no polys). Then I installed Houdini and I could import those files to the geometry view of Houdini. Now, I don't understand how we can get similar output video as shown in the repo (beautiful with plates and knife).

  • Correct implementation of C

    Correct implementation of C

    This computation of C as following is not consistent with this line. https://github.com/yuanming-hu/taichi_mpm/blob/3bb90fbe4c901aafc048dbb2d8d8aa388226d011/mls-mpm88-explained.cpp#L171 The above line is lack of one 1/dx according to the last page of this slide. 图片

  • Installation Issues

    Installation Issues

    Hi, I am trying to run the examples for mpm but I keep getting this error: Traceback (most recent call last): File "/taichi_mpm-master/scripts/mls-cpic/banana.py", line 8, in mpm = tc.dynamics.MPM( AttributeError: module 'taichi' has no attribute 'dynamics' I installed taichi using 'pip install'. Then I got the above error. I tried uninstalling and re-installing using the legacy link given but the install file keeps giving compile error. I have seen other people post about this as well under the "issues" section on github but no response. I hope you can help us out.

DeepLM: Large-scale Nonlinear Least Squares on Deep Learning Frameworks using Stochastic Domain Decomposition (CVPR 2021)
DeepLM: Large-scale Nonlinear Least Squares on Deep Learning Frameworks using Stochastic Domain Decomposition (CVPR 2021)

DeepLM DeepLM: Large-scale Nonlinear Least Squares on Deep Learning Frameworks using Stochastic Domain Decomposition (CVPR 2021) Run Please install th

May 10, 2022
T-LOAM: Truncated Least Squares Lidar-only Odometry and Mapping in Real-Time
T-LOAM: Truncated Least Squares Lidar-only Odometry and Mapping in Real-Time

T-LOAM: Truncated Least Squares Lidar-only Odometry and Mapping in Real-Time The first Lidar-only odometry framework with high performance based on tr

May 5, 2022
Official PyTorch implementation of the paper "Deep Constrained Least Squares for Blind Image Super-Resolution", CVPR 2022.
Official PyTorch implementation of the paper

Deep Constrained Least Squares for Blind Image Super-Resolution [Paper] This is the official implementation of 'Deep Constrained Least Squares for Bli

May 21, 2022
Simple Linear 2nd ODE Solver GUI - A 2nd constant coefficient linear ODE solver with simple GUI using euler's method
Simple Linear 2nd ODE Solver GUI - A 2nd constant coefficient linear ODE solver with simple GUI using euler's method

Simple_Linear_2nd_ODE_Solver_GUI Description It is a 2nd constant coefficient li

Feb 5, 2022
Sudoku solver - A sudoku solver with python

sudoku_solver A sudoku solver What is Sudoku? Sudoku (Japanese: 数独, romanized: s

Jan 6, 2022
A customisable game where you have to quickly click on black tiles in order of appearance while avoiding clicking on white squares.

W.I.P-Aim-Memory-Game A customisable game where you have to quickly click on black tiles in order of appearance while avoiding clicking on white squar

Dec 8, 2021
Creating a Linear Program Solver by Implementing the Simplex Method in Python with NumPy
Creating a Linear Program Solver by Implementing the Simplex Method in Python with NumPy

Creating a Linear Program Solver by Implementing the Simplex Method in Python with NumPy Simplex Algorithm is a popular algorithm for linear programmi

Jul 13, 2021
Feb 18, 2022
Rasterize with the least efforts for researchers.

utils3d Rasterize and do image-based 3D transforms with the least efforts for researchers. Based on numpy and OpenGL. It could be helpful when you wan

Apr 24, 2022
Style-based Point Generator with Adversarial Rendering for Point Cloud Completion (CVPR 2021)
Style-based Point Generator with Adversarial Rendering for Point Cloud Completion (CVPR 2021)

Style-based Point Generator with Adversarial Rendering for Point Cloud Completion (CVPR 2021) An efficient PyTorch library for Point Cloud Completion.

May 10, 2022
Implementation of the "PSTNet: Point Spatio-Temporal Convolution on Point Cloud Sequences" paper.
Implementation of the

PSTNet: Point Spatio-Temporal Convolution on Point Cloud Sequences Introduction Point cloud sequences are irregular and unordered in the spatial dimen

May 15, 2022
Code for "PV-RAFT: Point-Voxel Correlation Fields for Scene Flow Estimation of Point Clouds", CVPR 2021
Code for

PV-RAFT This repository contains the PyTorch implementation for paper "PV-RAFT: Point-Voxel Correlation Fields for Scene Flow Estimation of Point Clou

Apr 28, 2022
Implementation of the "Point 4D Transformer Networks for Spatio-Temporal Modeling in Point Cloud Videos" paper.
Implementation of the

Point 4D Transformer Networks for Spatio-Temporal Modeling in Point Cloud Videos Introduction Point cloud videos exhibit irregularities and lack of or

May 20, 2022
Point Cloud Denoising input segmentation output raw point-cloud valid/clear fog rain de-noised Abstract Lidar sensors are frequently used in environme
Point Cloud Denoising input segmentation output raw point-cloud valid/clear fog rain de-noised Abstract Lidar sensors are frequently used in environme

Point Cloud Denoising input segmentation output raw point-cloud valid/clear fog rain de-noised Abstract Lidar sensors are frequently used in environme

May 17, 2022
Synthetic LiDAR sequential point cloud dataset with point-wise annotations
Synthetic LiDAR sequential point cloud dataset with point-wise annotations

SynLiDAR dataset: Learning From Synthetic LiDAR Sequential Point Cloud This is official repository of the SynLiDAR dataset. For technical details, ple

May 12, 2022
[ICCV 2021 Oral] SnowflakeNet: Point Cloud Completion by Snowflake Point Deconvolution with Skip-Transformer
[ICCV 2021 Oral] SnowflakeNet: Point Cloud Completion by Snowflake Point Deconvolution with Skip-Transformer

This repository contains the source code for the paper SnowflakeNet: Point Cloud Completion by Snowflake Point Deconvolution with Skip-Transformer (ICCV 2021 Oral). The project page is here.

May 19, 2022
Unofficial implementation of Point-Unet: A Context-Aware Point-Based Neural Network for Volumetric Segmentation

Point-Unet This is an unofficial implementation of the MICCAI 2021 paper Point-Unet: A Context-Aware Point-Based Neural Network for Volumetric Segment

May 2, 2022
Point-NeRF: Point-based Neural Radiance Fields
Point-NeRF: Point-based Neural Radiance Fields

Point-NeRF: Point-based Neural Radiance Fields Project Sites | Paper | Primary c

May 18, 2022
DeFMO: Deblurring and Shape Recovery of Fast Moving Objects (CVPR 2021)
 DeFMO: Deblurring and Shape Recovery of Fast Moving Objects (CVPR 2021)

Evaluation, Training, Demo, and Inference of DeFMO DeFMO: Deblurring and Shape Recovery of Fast Moving Objects (CVPR 2021) Denys Rozumnyi, Martin R. O

May 17, 2022