Add dynamic levelset into 3D MPM simulator

include/taichi/dynamics/mpm2d/mpm.cpp

@@ -188,7 +188,6 @@ void MPM::substep() {
     estimate_volume();
     grid.backup_velocity();
     grid.apply_external_force(gravity);
-    //Deleted: grid.apply_boundary_conditions(levelset);
     apply_deformation_force();
     grid.normalize_acceleration();
     grid.apply_boundary_conditions(levelset, t_int_increment * base_delta_t, t);
@@ -230,9 +229,11 @@ void MPM::compute_material_levelset() {
 }
 
 void MPM::particle_collision_resolution() {
-    for (auto &p : particles) {
-        if (p->state == MPMParticle::UPDATING)
-            p->resolve_collision(levelset, t);
+    if (levelset.levelset0) {
+        for (auto &p : particles) {
+            if (p->state == MPMParticle::UPDATING)
+                p->resolve_collision(levelset, t);
+        }
     }
 }
 

include/taichi/dynamics/mpm2d/mpm_grid.cpp

@@ -14,7 +14,7 @@ TC_NAMESPACE_BEGIN
 long long MPMParticle::instance_count = 0;
 
 void Grid::apply_boundary_conditions(const DynamicLevelSet2D & levelset, real delta_t, real t) {
-    if (levelset.levelset0->get_width() > 0) {
+    if (levelset.levelset0) {
         for (auto &ind : boundary_normal.get_region()) {
             Vector2 pos = Vector2(ind.get_pos());
             Vector2 v = velocity[ind] + force_or_acc[ind] * delta_t - levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
@@ -28,7 +28,10 @@ void Grid::apply_boundary_conditions(const DynamicLevelSet2D & levelset, real de
                     v = Vector2(0.0f);
                 }
                 else {
-                    Vector2 t = Vector2(-n.y, n.x);
+                    Vector2 t = v - n * glm::dot(v, n);
+                    if (length(t) > 1e-6f) {
+                        t = normalize(t);
+                    }
                     real friction = -clamp(glm::dot(t, v), -mu * pressure, mu * pressure);
                     v = v + n * pressure + t * friction;
                 }

include/taichi/dynamics/simulation3d.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.
@@ -12,13 +13,15 @@
 #include <taichi/common/meta.h>
 #include <taichi/visualization/particle_visualization.h>
 #include <vector>
+#include <taichi/math/dynamic_levelset_3d.h>
 
 TC_NAMESPACE_BEGIN
 
 class Simulation3D : public Unit {
 protected:
     real current_t = 0.0f;
     int num_threads;
+    DynamicLevelSet3D levelset;
 public:
     Simulation3D() {}
     virtual real get_current_time() const {
@@ -34,6 +37,11 @@ class Simulation3D : public Unit {
         error("no impl");
         return std::vector<RenderParticle>();
     }
+
+    virtual void set_levelset(const DynamicLevelSet3D &levelset) {
+        this->levelset = levelset;
+    }
+
     virtual void update(const Config &config) {}
 };
 

include/taichi/math/dynamic_levelset_3d.cpp

@@ -32,12 +32,20 @@ Vector3 DynamicLevelSet3D::get_spatial_gradient(const Vector3 &pos, real t) cons
 }
 
 real DynamicLevelSet3D::get_temporal_derivative(const Vector3 &pos, real t) const {
+    assert_info(levelset0->inside(pos), "Temporal("
+                                        + std::to_string(pos.x) + ", "
+                                        + std::to_string(pos.y) + ", "
+                                        + std::to_string(pos.z) + ")");
     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 DynamicLevelSet3D::sample(const Vector3 &pos, real t) const {
+    assert_info(levelset0->inside(pos), "Sample("
+                             + std::to_string(pos.x) + ", "
+                             + std::to_string(pos.y) + ", "
+                             + std::to_string(pos.z) + ")");
     real l1 = levelset0->get(pos);
     real l2 = levelset1->get(pos);
     return lerp((t - t0) / (t1 - t0), l1, l2);
@@ -46,7 +54,7 @@ real DynamicLevelSet3D::get_temporal_derivative(const Vector3 &pos, real t) cons
 Array3D<real> DynamicLevelSet3D::rasterize(int width, int height, int depth, real t) {
     Array3D<real> r0 = levelset0->rasterize(width, height, depth);
     Array3D<real> r1 = levelset1->rasterize(width, height, depth);
-    Array3D<real> out(width, height, width);
+    Array3D<real> out(width, height, depth);
     for (auto &ind : Region3D(0, width, 0, height, 0, depth, Vector3(0.5f, 0.5f, 0.5f))) {
         out[ind] = lerp((t - t0) / (t1 - t0), r0[ind], r1[ind]);
         if (std::isnan(out[ind])) {

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

@@ -5,7 +5,7 @@
 
 if __name__ == '__main__':
     resolution = tuple([256] * 2)
-    simulator = create_mpm_simulator(resolution, 10, frame_dt=0.01)
+    simulator = create_mpm_simulator(resolution, 1, frame_dt=0.05)
 
     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))
 

python/examples/simulation/3d/mpm_snow_column_collapse_3d.py

@@ -0,0 +1,103 @@
+import math
+
+from taichi.dynamics.mpm import MPM3
+from taichi.core import tc_core
+from taichi.misc.util import Vector
+from taichi.visual import *
+from taichi.visual.post_process import *
+from taichi.visual.texture import Texture
+import taichi as tc
+
+gi_render = True
+step_number = 400
+# step_number = 1
+# total_frames = 1
+grid_downsample = 8
+output_downsample = 1
+render_epoch = 20
+
+
+def create_mpm_sand_block(fn):
+    particles = tc_core.RenderParticles()
+    assert particles.read(fn)
+    downsample = grid_downsample
+    tex = Texture.from_render_particles((255 / downsample, 255 / downsample, 255 / downsample), particles) * 5
+    # mesh_transform = tc_core.Matrix4(1.0).scale_s(0.5).translate(Vector(0.5, 0.5, 0.5))
+    # transform = tc_core.Matrix4(1.0).scale_s(2).scale(Vector(2.0, 0.5, 1.0)).translate(Vector(-2, -0.99, -1))
+    mesh_transform = tc_core.Matrix4(1.0).translate(Vector(0, 0.01, 0))
+    transform = tc_core.Matrix4(1.0).scale_s(2).translate(Vector(-1, -1, -1))
+    vol = VolumeMaterial('sdf_voxel', scattering=5, absorption=0, tex=tex,
+                         resolution=(255 / downsample, 255 / downsample, 255 / downsample),
+                         transform_ptr=transform.get_ptr_string())
+    material = SurfaceMaterial('plain_interface')
+    material.set_internal_material(vol)
+    return Mesh('cube', material=material, transform=transform * mesh_transform)
+
+
+def create_sand_scene(frame, d, t):
+    downsample = output_downsample
+    width, height = 540 / downsample, 540 / downsample
+    camera = Camera('thinlens', width=width, height=height, fov=75,
+                    origin=(0, 1, 4), look_at=(0.0, -0.9, 0.0), up=(0, 1, 0), aperture=0.01)
+
+    scene = Scene()
+    with scene:
+        scene.set_camera(camera)
+        rep = Texture.create_taichi_wallpaper(10, rotation=0, scale=0.95) * Texture('const', value=(0.7, 0.5, 0.5))
+        material = SurfaceMaterial('pbr', diffuse_map=rep)
+        scene.add_mesh(Mesh('holder', material=material, translate=(0, -1, -6), scale=2))
+
+        mesh = Mesh('plane', SurfaceMaterial('emissive', color=(1, 1, 1)),
+                    translate=(1.0, 1.0, -1), scale=(0.1, 0.1, 0.1), rotation=(180, 0, 0))
+        scene.add_mesh(mesh)
+
+        material = tc.SurfaceMaterial('microfacet', color=(1, 1, 0.5), roughness=(0.1, 0, 0, 0), f0=1)
+        sphere = tc.Mesh('sphere', material,
+                         translate=((t+0.05) * 0.5 - 0.35, -0.61, 0), scale=0.1, rotation=(0, 0, 0))
+        scene.add_mesh(sphere)
+
+        # Change this line to your particle output path pls.
+        # fn = r'../sand-sim/particles%05d.bin' % frame
+        fn = d + r'/particles%05d.bin' % frame
+        mesh = create_mpm_sand_block(fn)
+        scene.add_mesh(mesh)
+
+    return scene
+
+
+def render_sand_frame(frame, d, t):
+    renderer = Renderer(output_dir='volumetric', overwrite=True, frame=frame)
+    renderer.initialize(preset='pt', scene=create_sand_scene(frame, d, t), sampler='prand')
+    renderer.set_post_processor(LDRDisplay(exposure=0.6, bloom_radius=0.0, bloom_threshold=1.0))
+    renderer.render(render_epoch)
+
+
+if __name__ == '__main__':
+    downsample = grid_downsample
+    resolution = (255 / downsample, 255 / downsample, 255 / downsample)
+    tex = Texture('ring', outer=0.15) * 4
+    tex = Texture('bound', tex=tex, axis=2, bounds=(0.0, 0.4), outside_val=(0, 0, 0))
+    tex = Texture('rotate', tex=tex, rotate_axis=0, rotate_times=1)
+    mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), initial_velocity=(0, 0, 0), delta_t=0.00005, num_threads=8,
+               density_tex=tex.id)
+
+
+    # Dynamic Levelset
+    def levelset_generator(t):
+        levelset = mpm.create_levelset()
+        levelset.add_sphere(Vector(0.325 + 0.25 * (t+0.05), 0.2, 0.5), 0.1, False)
+        # levelset.add_sphere(Vector(0.5, 0.2, 0.5), t, False)
+        return levelset
+    mpm.set_levelset(levelset_generator, True)
+
+    t = 0
+    for i in range(step_number):
+        print 'process(%d/%d)' % (i, step_number)
+        mpm.step(0.01)
+        t += 0.01
+        if gi_render:
+            d = mpm.get_directory()
+            if i % 10 == 0:
+                render_sand_frame(i, d, t)
+                pass
+    mpm.make_video()

python/taichi/dynamics/levelset_3d.py

@@ -0,0 +1,31 @@
+import taichi as tc
+from taichi.misc.util import *
+
+
+class LevelSet3D:
+    def __init__(self, res, offset=None):
+        if offset is None:
+            offset = Vector(0, 0, 0)
+        self.delta_x = 1.0 / min(res)
+        self.res = (res[0] + 1, res[1] + 1, res[2] + 1)
+        self.levelset = tc.core.LevelSet3D(int(res[0]), int(res[1]), int(res[2]), offset)
+
+    def add_sphere(self, center, radius, inside_out=False):
+        if type(center) != tc.core.Vector3:
+            center = Vector(center[0], center[1], center[2])
+        self.levelset.add_sphere(Vector(center.x / self.delta_x, center.y / self.delta_x, center.z / self.delta_x),
+                                 radius / self.delta_x, inside_out)
+
+    # def get(self, x, y=None):
+    #     if y is None:
+    #         y = x.y
+    #         x = x.x
+    #     return self.levelset.sample(x / self.delta_x, y / self.delta_x)
+    #
+    # def get_normalized_gradient(self, p):
+    #     p.x /= self.delta_x
+    #     p.y /= self.delta_x
+    #     return self.levelset.get_normalized_gradient(p)
+
+    def set_friction(self, f):
+        self.levelset.friction = f

python/taichi/dynamics/mpm.py

@@ -1,6 +1,7 @@
 import time
 
 from taichi.core import tc_core
+from levelset_3d import LevelSet3D
 from taichi.misc.util import *
 from taichi.tools.video import VideoManager
 from taichi.visual.camera import Camera
@@ -27,11 +28,29 @@ def __init__(self, **kwargs):
                                                   light_direction=(1, 1, 0))
         self.resolution = kwargs['resolution']
         self.frame = 0
+        self.levelset_generator = None
+
+    def update_levelset(self, t0, t1):
+        levelset = tc.core.DynamicLevelSet3D()
+        levelset.initialize(t0, t1, self.levelset_generator(t0).levelset, self.levelset_generator(t1).levelset)
+        self.c.set_levelset(levelset)
+
+    def set_levelset(self, levelset, is_dynamic_levelset = False):
+        if is_dynamic_levelset:
+            self.levelset_generator = levelset
+        else:
+            def levelset_generator(_):
+                return levelset
+            self.levelset_generator = levelset_generator
+
+    def get_current_time(self):
+        return self.c.get_current_time()
 
     def step(self, step_t):
         t = self.c.get_current_time()
         print 'Simulation time:', t
         T = time.time()
+        self.update_levelset(t, t + step_t)
         self.c.step(step_t)
         print 'Time:', time.time() - T
         image_buffer = tc_core.Array2DVector3(self.video_manager.width, self.video_manager.height, Vector(0, 0, 0.0))
@@ -54,3 +73,6 @@ def get_directory(self):
 
     def make_video(self):
         self.video_manager.make_video()
+
+    def create_levelset(self):
+        return LevelSet3D(self.resolution, Vector(0.0, 0.0, 0.0))

python/taichi/two_d/levelset_2d.py

@@ -1,10 +1,11 @@
 import taichi as tc
 from taichi.misc.util import *
 
+
 class LevelSet2D:
     def __init__(self, width, height, delta_x, offset=None):
         if offset is None:
-            offset = Vector(0.5, 0.5)
+            offset = Vector(0, 0)
         self.delta_x = delta_x
         self.levelset = tc.core.LevelSet2D(width, height, offset)
         self.width = width
@@ -20,16 +21,16 @@ def add_sphere(self, center, radius, inside_out=False):
     def add_polygon(self, polygon, inside_out):
         self.levelset.add_polygon(make_polygon(polygon, 1.0 / self.delta_x), inside_out)
 
-    def get(self, x, y=None):
-        if y is None:
-            y = x.y
-            x = x.x
-        return self.levelset.sample(x / self.delta_x, y / self.delta_x)
-
-    def get_normalized_gradient(self, p):
-        p.x /= self.delta_x
-        p.y /= self.delta_x
-        return self.levelset.get_normalized_gradient(p)
+    # def get(self, x, y=None):
+    #     if y is None:
+    #         y = x.y
+    #         x = x.x
+    #     return self.levelset.sample(x / self.delta_x, y / self.delta_x)
+    #
+    # def get_normalized_gradient(self, p):
+    #     p.x /= self.delta_x
+    #     p.y /= self.delta_x
+    #     return self.levelset.get_normalized_gradient(p)
 
     def get_image(self, width, height, color_255):
         if self.cache_image is None:

python/taichi/two_d/simulation_window.py

@@ -38,6 +38,7 @@ def __init__(self, max_side, simulator, color_scheme, levelset_supersampling=2,
         self.levelset_supersampling = levelset_supersampling
         self.show_grid = show_grid
         self.quit_pressed = False
+        self.timer = 0
         pyglet.clock.schedule_interval(self.update, 1 / 120.0)
         pyglet.app.run()
 
@@ -56,16 +57,19 @@ def clear_frames(self):
         os.rmdir(self.frame_output_path)
 
     def update(self, _):
+        self.timer += 1
+        print self.timer
+        if self.timer % 5 == 0:
+            self.make_video()
         if not self.quit_pressed and not self.simulator.ended():
             self.simulator.step()
         else:
-            self.make_video()
             exit(0)
         self.redraw()
         self.save_frame()
 
     def make_video(self):
-        command = "ffmpeg -framerate 24 -i " + self.frame_output_path + \
+        command = "ffmpeg -y -framerate 24 -i " + self.frame_output_path + \
                   "/frame%d.png -s:v " + str(self.width) + 'x' + str(self.height) + \
                   " -c:v libx264 -profile:v high -crf 20 -pix_fmt yuv420p " + VIDEO_OUTPUT_ROOT + \
                   "/" + self.video_filename

src/python/export_dynamics.cpp

@@ -72,6 +72,7 @@ void export_dynamics(py::module &m) {
         .def("step", &SIM::step) \
         .def("get_current_time", &SIM::get_current_time) \
         .def("get_render_particles", &SIM::get_render_particles) \
+        .def("set_levelset", &SIM::set_levelset) \
         ;
     EXPORT_SIMULATOR_3D(Simulation3D);