Add some tests and fix some bugs for 3D sand simulation

python/examples/rendering/test_volumetric.py

@@ -82,5 +82,5 @@ def render_snow_frame(frame):
 
 if __name__ == '__main__':
     total_frames = 550
-    for i in range(0, total_frames, 5):
+    for i in range(0, 1, 5):
         render_snow_frame(i)

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

@@ -4,7 +4,7 @@
 
 if __name__ == '__main__':
     downsample = 2
-    resolution = (255 / downsample, 255 / downsample, 255 / downsample)
+    resolution = (511 / downsample, 127 / downsample, 255 / downsample)
     tex = Texture('image', filename=tc.get_asset_path('/textures/taichi_words.png')) * 8
     tex = Texture('bound', tex=tex, axis=2, bounds=(0.475, 0.525), outside_val=(0, 0, 0))
     mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), initial_velocity=(0, 0, 0), delta_t=0.001, num_threads=8,

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

@@ -0,0 +1,76 @@
+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
+
+step_number = 6000
+# step_number = 1
+# total_frames = 1
+grid_downsample = 1
+output_downsample = 1
+render_epoch = 20
+
+
+def create_mpm_snow_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))
+    vol = VolumeMaterial('sdf_voxel', scattering=5, absorption=0, tex=tex, resolution=(255, 255, 255),
+                         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_snow_scene(frame, d):
+    downsample = output_downsample
+    width, height = 540 / downsample, 540 / downsample
+    camera = Camera('thinlens', width=width, height=height, fov=35,
+                    origin=(0, 1, 4), look_at=(0.0, -0.9, -0.0), up=(0, 1, 0), aperture=0.08)
+
+    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)
+
+        # Change this line to your particle output path pls.
+        # fn = r'../snow-sim/particles%05d.bin' % frame
+        fn = d + r'/particles%05d.bin' % frame
+        mesh = create_mpm_snow_block(fn)
+        scene.add_mesh(mesh)
+
+    return scene
+
+
+def render_snow_frame(frame, d):
+    renderer = Renderer(output_dir='volumetric', overwrite=True, frame=frame)
+    renderer.initialize(preset='pt', scene=create_snow_scene(frame, d), 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.038) * 8
+    tex = Texture('bound', tex=tex, axis=2, bounds=(0.0, 0.35), 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.0005, num_threads=8,
+               density_tex=tex.id)
+    for i in range(step_number):
+        print 'process(%d/%d)' % (i, step_number)
+        mpm.step(0.01)
+        d = mpm.get_directory()
+        if i % 20 == 0:
+            render_snow_frame(i, d)

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

@@ -0,0 +1,80 @@
+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
+
+step_number = 6000
+# step_number = 1
+# total_frames = 1
+grid_downsample = 16
+output_downsample = 1
+render_epoch = 20
+grid_x = 255 / grid_downsample
+grid_y = 255 / grid_downsample
+grid_z = 255 / grid_downsample
+
+
+def create_mpm_snow_block(fn):
+    particles = tc_core.RenderParticles()
+    assert particles.read(fn)
+    tex = Texture.from_render_particles((grid_x, grid_y, grid_z), particles) * 5
+    mesh_transform = tc_core.Matrix4(1.0).scale_s(0.3).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))
+    vol = VolumeMaterial('sdf_voxel', scattering=5, absorption=0, tex=tex, resolution=(grid_x, grid_y, grid_z),
+                         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_snow_scene(frame, d):
+    downsample = output_downsample
+    width, height = 540 / downsample, 540 / downsample
+    camera = Camera('thinlens', width=width, height=height, fov=90,
+                    origin=(0, 1, 4), look_at=(0.0, -0.9, -0.0), up=(0, 1, 0), aperture=0.08)
+
+    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)
+
+        # Change this line to your particle output path pls.
+        # fn = r'../snow-sim/particles%05d.bin' % frame
+        fn = d + r'/particles%05d.bin' % frame
+        mesh = create_mpm_snow_block(fn)
+        scene.add_mesh(mesh)
+
+    return scene
+
+
+def render_snow_frame(frame, d):
+    renderer = Renderer(output_dir='volumetric', overwrite=True, frame=frame)
+    renderer.initialize(preset='pt', scene=create_snow_scene(frame, d), 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 = (grid_x, grid_y, grid_z)
+    # tex = Texture('ring', outer=0.038) * 8
+    tex = Texture('rect', bounds=(1, 1, 1)) * 8
+    tex = Texture('bound', tex=tex, axis=2, bounds=(0.0, 0.5), 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.0002, num_threads=8,
+               density_tex=tex.id)
+    for i in range(step_number):
+        print 'process(%d/%d)' % (i, step_number)
+        # mpm.step(0.01)
+        d = mpm.get_directory()
+        d = r'/Users/squarefk/repos/taichi_outputs/task-2017-03-30-13-51-25-r09595'
+        if i % 20 == 0:
+            render_snow_frame(i, d)

python/taichi/dynamics/mpm.py

@@ -46,5 +46,8 @@ def step(self, step_t):
         self.video_manager.write_frame(img)
         self.frame += 1
 
+    def get_directory(self):
+        return self.directory
+
     def make_video(self):
         self.video_manager.make_video()

src/simulation3d/mpm/mpm3.cpp

@@ -99,17 +99,8 @@ struct EPParticle3 : public MPM3D::Particle {
 };
 
 struct DPParticle3 : public MPM3D::Particle {
-//    particle.mu_0 = 1e6
-//    particle.lambda_0 = 2e5
-//    particle.mu_0 = 10000000
-//    particle.lambda_0 = 10000000
-//    particle.h_0 = 45
-//    particle.h_1 = 9
-//    particle.h_2 = 0.2
-//    particle.h_3 = 10
-//    particle.alpha = 1
     real h_0 = 35, h_1 = 9, h_2 = 0.2, h_3 = 10;
-    real lambda_0 = 200000, mu_0 = 1000000;
+    real lambda_0 = 204057, mu_0 = 136038;
     real alpha = 1;
     real q = 0.0f;
     real phi_f;

src/texture/texture_op.cpp

@@ -130,19 +130,33 @@ TC_IMPLEMENTATION(Texture, RepeatedTexture, "repeat");
 class RotatedTexture : public Texture {
 protected:
     std::shared_ptr<Texture> tex;
-    int times;
+    int rotate_times;
+    int rotate_axis;
 public:
     void initialize(const Config &config) override {
         Texture::initialize(config);
         tex = AssetManager::get_asset<Texture>(config.get_int("tex"));
-        times = config.get_int("times");
+        rotate_times = config.get_int("rotate_times");
+        rotate_axis = config.get_int("rotate_axis");
     }
 
     virtual Vector4 sample(const Vector3 &coord_) const override {
         auto coord = coord_;
-        for (int i = 0; i < times; i++) {
-            coord = Vector3(-coord.y, coord.x, coord.z);
+        coord = coord * 2.f - Vector3(1.f, 1.f, 1.f);
+        for (int i = 0; i < rotate_times; i++) {
+            switch (rotate_axis) {
+                case 0:
+                    coord = Vector3(coord.x, -coord.z, coord.y);
+                    break;
+                case 1:
+                    coord = Vector3(coord.z, coord.y, -coord.x);
+                    break;
+                case 2:
+                    coord = Vector3(-coord.y, coord.x, coord.z);
+                    break;
+            }
         }
+        coord = (coord + Vector3(1.f, 1.f, 1.f)) / 2.f;
         return tex->sample(coord);
     }
 };