Replace static levelset with dynamic levelset for 2D mpm simulator

include/taichi/dynamics/mpm2d/mpm.cpp

@@ -2,6 +2,7 @@
     Taichi - Physically based Computer Graphics Library
 
     Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
+                  2017 Yu Fang <squarefk@gmail.com>
 
     All rights reserved. Use of this source code is governed by
     the MIT license as written in the LICENSE file.
@@ -130,7 +131,7 @@ void MPM::substep() {
     //Deleted: grid.apply_boundary_conditions(levelset);
     apply_deformation_force();
     grid.normalize_acceleration();
-    grid.apply_boundary_conditions(levelset, t_int_increment * base_delta_t);
+    grid.apply_boundary_conditions(levelset, t_int_increment * base_delta_t, t);
     resample();
     for (auto &p : particles) {
         if (p->state == MPMParticle::UPDATING) {
@@ -162,7 +163,7 @@ void MPM::compute_material_levelset() {
     }
     for (auto &ind : material_levelset.get_region()) {
         if (material_levelset[ind] < 0.5f) {
-            if (levelset.sample(ind.get_pos()) < 0)
+            if (levelset.sample(ind.get_pos(), t) < 0)
                 material_levelset[ind] = -0.5f;
         }
     }
@@ -171,7 +172,7 @@ void MPM::compute_material_levelset() {
 void MPM::particle_collision_resolution() {
     for (auto &p : particles) {
         if (p->state == MPMParticle::UPDATING)
-            p->resolve_collision(levelset);
+            p->resolve_collision(levelset, t);
     }
 }
 

include/taichi/dynamics/mpm2d/mpm.h

@@ -2,6 +2,7 @@
     Taichi - Physically based Computer Graphics Library
 
     Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
+                  2017 Yu Fang <squarefk@gmail.com>
 
     All rights reserved. Use of this source code is governed by
     the MIT license as written in the LICENSE file.
@@ -13,6 +14,7 @@
 #include <vector>
 #include "mpm_grid.h"
 #include <taichi/math/levelset_2d.h>
+#include <taichi/math/dynamic_levelset_2d.h>
 #include <taichi/visual/texture.h>
 #include <taichi/visualization/image_buffer.h>
 
@@ -42,7 +44,7 @@ class MPM {
     Vector2 gravity;
     bool apic;
 
-    LevelSet2D levelset;
+    DynamicLevelSet2D levelset;
     LevelSet2D material_levelset;
 
     real last_sort;
@@ -97,7 +99,7 @@ class MPM {
 
     real get_current_time();
 
-    void set_levelset(const LevelSet2D &levelset) {
+    void set_levelset(const DynamicLevelSet2D &levelset) {
         this->levelset = levelset;
     }
 

include/taichi/dynamics/mpm2d/mpm_grid.cpp

@@ -13,15 +13,17 @@ TC_NAMESPACE_BEGIN
 
 long long MPMParticle::instance_count = 0;
 
-void Grid::apply_boundary_conditions(const LevelSet2D & levelset, real delta_t) {
-    if (levelset.get_width() > 0) {
+void Grid::apply_boundary_conditions(const DynamicLevelSet2D & levelset, real delta_t, real t) {
+    if (levelset.levelset0->get_width() > 0) {
         for (auto &ind : boundary_normal.get_region()) {
-            Vector2 v = velocity[ind] + force_or_acc[ind] * delta_t, n = levelset.get_normalized_gradient(Vector2(ind.i + 0.5f, ind.j + 0.5f));
-            real phi = levelset[ind];
+            Vector2 pos = Vector2(ind.i + 0.5f, ind.j + 0.5f);
+            Vector2 v = velocity[ind] + force_or_acc[ind] * delta_t - levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
+            Vector2 n = levelset.get_spatial_gradient(pos, t);
+            real phi = levelset.sample(pos, t);
             if (phi > 1) continue;
             else if (phi > 0) { // 0~1
                 real pressure = std::max(-glm::dot(v, n), 0.0f);
-                real mu = levelset.friction;
+                real mu = levelset.levelset0->friction;
                 if (mu < 0) { // sticky
                     v = Vector2(0.0f);
                 }
@@ -34,6 +36,7 @@ void Grid::apply_boundary_conditions(const LevelSet2D & levelset, real delta_t)
             else if (phi <= 0) {
                 v = Vector2(0.0f);
             }
+            v += levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
             force_or_acc[ind] = (v - velocity[ind]) / delta_t;
         }
     }

include/taichi/dynamics/mpm2d/mpm_grid.h

@@ -2,6 +2,7 @@
     Taichi - Physically based Computer Graphics Library
 
     Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
+                  2017 Yu Fang <squarefk@gmail.com>
 
     All rights reserved. Use of this source code is governed by
     the MIT license as written in the LICENSE file.
@@ -16,6 +17,7 @@
 #include <taichi/math/array_2d.h>
 #include <taichi/math/array_1d.h>
 #include <taichi/math/levelset_2d.h>
+#include <taichi/math/dynamic_levelset_2d.h>
 #include "mpm_particle.h"
 
 TC_NAMESPACE_BEGIN
@@ -135,7 +137,7 @@ class Grid {
         }
     }
 
-    void apply_boundary_conditions(const LevelSet2D &levelset, real delta_t);
+    void apply_boundary_conditions(const DynamicLevelSet2D &levelset, real delta_t, real t);
 
     void check_velocity() {
         for (int i = 0; i < res[0]; i++) {

include/taichi/dynamics/mpm2d/mpm_particle.h

@@ -2,6 +2,7 @@
     Taichi - Physically based Computer Graphics Library
 
     Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
+                  2017 Yu Fang <squarefk@gmail.com>
 
     All rights reserved. Use of this source code is governed by
     the MIT license as written in the LICENSE file.
@@ -14,6 +15,7 @@
 #include "mpm_utils.h"
 #include <taichi/common/util.h>
 #include <taichi/math/levelset_2d.h>
+#include <taichi/math/dynamic_levelset_2d.h>
 
 TC_NAMESPACE_BEGIN
 
@@ -104,10 +106,10 @@ struct MPMParticle {
 
     virtual void plasticity() {};
 
-    virtual void resolve_collision(const LevelSet2D &levelset) {
-        real phi = levelset.sample(pos.x, pos.y);
+    virtual void resolve_collision(const DynamicLevelSet2D &levelset, real t) {
+        real phi = levelset.sample(pos, t);
         if (phi < 0) {
-            pos -= levelset.get_normalized_gradient(pos) * phi;
+            pos -= levelset.get_spatial_gradient(pos, t) * phi;
         }
     }
 

include/taichi/math/dynamic_levelset_2d.cpp

@@ -18,18 +18,28 @@ void DynamicLevelSet2D::initialize(real _t0, real _t1, const LevelSet2D &_ls0, c
     levelset1 = std::make_shared<LevelSet2D>(_ls1);
 }
 
-Vector2 DynamicLevelSet2D::get_spatial_gradient(const Vector2 &pos, real t) {
+Vector2 DynamicLevelSet2D::get_spatial_gradient(const Vector2 &pos, real t) const {
     Vector2 gxy0 = levelset0->get_gradient(pos);
     Vector2 gxy1 = levelset1->get_gradient(pos);
     real gx = lerp((t - t0) / (t1 - t0), gxy0.x, gxy1.x);
     real gy = lerp((t - t0) / (t1 - t0), gxy0.y, gxy1.y);
-    return Vector2(gx, gy);
+    Vector2 gradient = Vector2(gx, gy);
+    if (length(gradient) < 1e-10f)
+        return Vector2(1, 0);
+    else
+        return normalize(gradient);
 }
 
-real DynamicLevelSet2D::get_temporal_gradient(const Vector2 &pos, real t) {
+real DynamicLevelSet2D::get_temporal_derivative(const Vector2 &pos, real t) const {
     real l1 = levelset0->get(pos);
-    real l2 = levelset1->get(pos);
-    return (l2 - l1) / (t1 - t0);
+    real l0 = levelset1->get(pos);
+    return (l1 - l0) / (t1 - t0);
+}
+
+real DynamicLevelSet2D::sample(const Vector2 &pos, real t) const {
+    real l1 = levelset0->sample(pos);
+    real l2 = levelset1->sample(pos);
+    return lerp((t - t0) / (t1 - t0), l1, l2);
 }
 
 Array2D<real> DynamicLevelSet2D::rasterize(int width, int height, real t) {

include/taichi/math/dynamic_levelset_2d.h

@@ -24,12 +24,14 @@ class DynamicLevelSet2D {
 
         void initialize(real _t0, real _t1, const LevelSet2D &_ls0, const LevelSet2D &_ls1);
 
-        Vector2 get_spatial_gradient(const Vector2 &pos, real t);
+        Vector2 get_spatial_gradient(const Vector2 &pos, real t) const;
 
-        real get_temporal_gradient(const Vector2 &pos, real t);
+        real get_temporal_derivative(const Vector2 &pos, real t) const;
+
+        real sample(const Vector2 &pos, real t) const;
 
         Array2D<real> rasterize(int width, int height, real t);
-    };
+};
 
 TC_NAMESPACE_END
 

include/taichi/math/dynamic_levelset_3d.cpp

@@ -18,19 +18,29 @@ void DynamicLevelSet3D::initialize(real _t0, real _t1, const LevelSet3D &_ls0, c
     levelset1 = std::make_shared<LevelSet3D>(_ls1);
 }
 
-Vector3 DynamicLevelSet3D::get_spatial_gradient(const Vector3 &pos, real t) {
+Vector3 DynamicLevelSet3D::get_spatial_gradient(const Vector3 &pos, real t) const {
     Vector3 gxyz0 = levelset0->get_gradient(pos);
     Vector3 gxyz1 = levelset1->get_gradient(pos);
     real gx = lerp((t - t0) / (t1 - t0), gxyz0.x, gxyz1.x);
     real gy = lerp((t - t0) / (t1 - t0), gxyz0.y, gxyz1.y);
     real gz = lerp((t - t0) / (t1 - t0), gxyz0.z, gxyz1.z);
-    return Vector3(gx, gy, gz);
+    Vector3 gradient = Vector3(gx, gy, gz);
+    if (length(gradient) < 1e-10f)
+        return Vector3(1, 0, 0);
+    else
+        return normalize(gradient);
 }
 
-real DynamicLevelSet3D::get_temporal_gradient(const Vector3 &pos, real t) {
+real DynamicLevelSet3D::get_temporal_derivative(const Vector3 &pos, real t) const {
+    real l0 = levelset0->get(pos);
+    real l1 = levelset1->get(pos);
+    return (l1 - l0) / (t1 - t0);
+}
+
+    real DynamicLevelSet3D::sample(const Vector3 &pos, real t) const {
     real l1 = levelset0->get(pos);
     real l2 = levelset1->get(pos);
-    return (l2 - l1) / (t1 - t0);
+    return lerp((t - t0) / (t1 - t0), l1, l2);
 }
 
 Array3D<real> DynamicLevelSet3D::rasterize(int width, int height, int depth, real t) {

include/taichi/math/dynamic_levelset_3d.h

@@ -24,9 +24,11 @@ class DynamicLevelSet3D {
 
     void initialize(real _t0, real _t1, const LevelSet3D &_ls0, const LevelSet3D &_ls1);
 
-    Vector3 get_spatial_gradient(const Vector3 &pos, real t);
+    Vector3 get_spatial_gradient(const Vector3 &pos, real t) const;
 
-    real get_temporal_gradient(const Vector3 &pos, real t);
+    real get_temporal_derivative(const Vector3 &pos, real t) const;
+
+    real sample(const Vector3 &pos, real t) const;
 
     Array3D<real> rasterize(int width, int height, int depth, real t);
 };

python/examples/simulation/2d/mpm_sand_stir.py

@@ -0,0 +1,28 @@
+import math
+import taichi as tc
+from taichi.two_d import *
+from taichi.misc.util import *
+
+if __name__ == '__main__':
+    resolution = tuple([256] * 2)
+    simulator = create_mpm_simulator(resolution, 10, frame_dt=0.01)
+
+    simulator.add_event(-1, lambda s: s.add_particles_polygon([(0.45, 0.15), (0.55, 0.15), (0.55, 0.8), (0.45, 0.8)], 'dp', h_0=20))
+
+    # Static Levelset
+    # levelset = simulator.create_levelset()
+    # levelset.add_polygon([(0.1, 0.1), (0.9, 0.1), (0.9, 0.9), (0.1, 0.9)], True)
+    # simulator.set_levelset(levelset)
+
+    # Dynamic Levelset
+    def levelset_generator(t):
+        levelset = simulator.create_levelset()
+        velocity = (0.5 + 0.25 * t) * 3.1415
+        print 0.25 * math.cos(t * velocity), 0.25 * math.sin(t * velocity)
+        levelset.add_sphere(Vector(0.5 + 0.25 * math.cos(t * velocity), 0.5 + 0.25 * math.sin(t * velocity)), 0.1, False)
+        levelset.add_sphere(Vector(0.5 + 0.25 * math.cos(t * velocity + math.pi), 0.5 + 0.25 * math.sin(t * velocity + math.pi)), 0.1, False)
+        levelset.add_sphere(Vector(0.5, 0.5), 0.45, True)
+        return levelset
+    simulator.set_levelset(levelset_generator, True)
+
+    window = SimulationWindow(512, simulator, color_schemes['sand'], levelset_supersampling=2, show_images=True)

python/taichi/two_d/mpm_simulator.py

@@ -146,8 +146,10 @@ def add_particles(self, particles):
 
     def get_levelset_images(self, width, height, color_scheme):
         images = []
-        images.append(self.levelset.get_image(width, height, color_scheme['boundary']))
-        material_levelset = self.get_material_levelset()
+        t = self.simulator.get_current_time()
+        levelset = self.levelset_generator(t)
+        images.append(levelset.get_image(width, height, color_scheme['boundary']))
+        material_levelset = self.simulator.get_material_levelset()
         images.append(array2d_to_image(material_levelset, width, height, color_scheme['material']))
         # if self.buffer_debug:
         #     images.append(array2d_to_image(material_levelset, width, height, color_scheme['material']))

python/taichi/two_d/simulator.py

@@ -7,6 +7,7 @@ def __init__(self, simulation_time, dt):
         self.events = []
         self.simulator = None
         self.simulation_time = simulation_time
+        self.levelset_generator = None
         self.dt = dt
         self.delta_x = 1
         self.particles = []
@@ -15,8 +16,14 @@ def add_event(self, t, func):
         self.events.append((t, func))
         self.events.sort()
 
+    def update_levelset(self, t0, t1):
+        levelset = tc.core.DynamicLevelSet2D()
+        levelset.initialize(t0, t1, self.levelset_generator(t0).levelset, self.levelset_generator(t1).levelset)
+        self.simulator.set_levelset(levelset)
+
     def step(self):
         t = self.simulator.get_current_time()
+        self.update_levelset(t, t + self.dt)
         while self.events and t > self.events[0][0]:
             self.events[0][1](self)
             self.events = self.events[1:]
@@ -39,13 +46,19 @@ def add_source(self, **kwargs):
     def ended(self):
         return self.simulator.get_current_time() >= self.simulation_time
 
-    def set_levelset(self, levelset):
-        self.simulator.set_levelset(levelset.levelset)
-        self.levelset = levelset
+    def set_levelset(self, levelset, is_dynamic_levelset = False):
+        if is_dynamic_levelset:
+            self.levelset_generator = levelset
+        else:
+            def levelset_generator(t):
+                return levelset
+            self.levelset_generator = levelset_generator
 
     def get_levelset_images(self, width, height, color_scheme):
         images = []
-        images.append(self.levelset.get_image(width, height, color_scheme['boundary']))
+        t = self.simulator.get_current_time()
+        levelset = self.levelset_generator(t)
+        images.append(levelset.get_image(width, height, color_scheme['boundary']))
         return images
 
     def maginify(self, val):

src/python/export_math.cpp

@@ -12,6 +12,7 @@
 #include <taichi/math/levelset_2d.h>
 #include <taichi/visualization/rgb.h>
 #include <taichi/math/array_op.h>
+#include <taichi/math/dynamic_levelset_2d.h>
 
 PYBIND11_MAKE_OPAQUE(std::vector<int>);
 PYBIND11_MAKE_OPAQUE(std::vector<taichi::real>);
@@ -111,6 +112,10 @@ void export_math(py::module &m) {
             .def("to_ndarray", &array2d_to_ndarray<LevelSet2D>)
             .def_readwrite("friction", &LevelSet2D::friction);
 
+    py::class_<DynamicLevelSet2D>(m, "DynamicLevelSet2D")
+            .def(py::init<>())
+            .def("initialize", &DynamicLevelSet2D::initialize);
+
     m.def("points_inside_polygon", points_inside_polygon);
     m.def("points_inside_sphere", points_inside_sphere);
     m.def("make_range", make_range);