Merge pull request #41 from squarefk/ampm

Throw balls scene. Thanks to @squarefk

include/taichi/math/levelset_3d.cpp

@@ -21,6 +21,21 @@ void LevelSet3D::add_sphere(Vector3 center, real radius, bool inside_out) {
     }
 }
 
+void LevelSet3D::add_plane(real a, real b, real c, real d) {
+    real coeff = 1.0f / sqrt(a * a + b * b + c * c);
+    for (auto &ind : get_region()) {
+        Vector3 sample = ind.get_pos();
+        real dist = (glm::dot(sample, Vector3(a, b, c)) + d) * coeff;
+        set(ind, std::min(Array3D::get(ind), dist));
+    }
+}
+
+void LevelSet3D::global_increase(real delta) {
+    for (auto &ind : get_region()) {
+        set(ind, Array3D::get(ind) + delta);
+    }
+}
+
 Vector3 LevelSet3D::get_gradient(const Vector3 &pos) const {
     assert_info(inside(pos), "LevelSet Gradient Query out of Bound! ("
                              + std::to_string(pos.x) + ", "

include/taichi/math/levelset_3d.h

@@ -45,6 +45,10 @@ class LevelSet3D : public Array3D<real> {
 
     void add_sphere(Vector3 center, real radius, bool inside_out = false);
 
+    void add_plane(real a, real b, real c, real d);
+
+    void global_increase(real delta);
+
     Vector3 get_gradient(const Vector3 &pos) const; // Note this is not normalized!
 
     Vector3 get_normalized_gradient(const Vector3 &pos) const;

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

@@ -0,0 +1,97 @@
+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 = False
+step_number = 400
+# step_number = 1
+# total_frames = 1
+grid_downsample = 3
+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, 3.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'../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)
+    levelset = MPM3.create_levelset(resolution)
+    levelset.add_plane(0, 1, 0, -1)
+    levelset.add_plane(1, 1, 0, -1.7)
+    levelset.global_increase(1)
+    tex = Texture('levelset3d', levelset=levelset, bounds=(0, 0.05 / levelset.get_delta_x())) * 8
+    tex = Texture('bound', tex=tex, axis=2, bounds=(0.3, 0.7), outside_val=(0, 0, 0))
+
+    tex_ball = Texture('sphere', center=(0.1, 0.5, 0.5), radius=0.05) * 8
+    tex += tex_ball
+
+    mpm = MPM3(resolution=resolution, gravity=(0, -100, 0), initial_velocity=(0, 0, 0), delta_t=0.001, num_threads=8,
+               density_tex=tex.id)
+    mpm.set_levelset(levelset, False)
+
+    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

@@ -9,13 +9,20 @@ def __init__(self, res, offset=None):
         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)
+        self.id = tc.core.register_levelset3d(self.levelset)
 
     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 add_plane(self, a, b, c, d):
+        self.levelset.add_plane(a, b, c, d / self.delta_x)
+
+    def global_increase(self, delta):
+        self.levelset.global_increase(delta / self.delta_x)
+
     # def get(self, x, y=None):
     #     if y is None:
     #         y = x.y
@@ -29,3 +36,6 @@ def add_sphere(self, center, radius, inside_out=False):
 
     def set_friction(self, f):
         self.levelset.friction = f
+
+    def get_delta_x(self):
+        return self.delta_x

python/taichi/dynamics/mpm.py

@@ -84,5 +84,9 @@ def make_video(self):
     def create_levelset(self):
         return LevelSet3D(self.resolution, Vector(0.0, 0.0, 0.0))
 
+    @staticmethod
+    def create_levelset(res):
+        return LevelSet3D(res, Vector(0.0, 0.0, 0.0))
+
     def test(self):
         return self.c.test()

python/taichi/visual/asset_manager.py

@@ -1,9 +1,11 @@
 import surface_material
 import volume_material
 import texture
+from taichi.dynamics import levelset_3d
 
 def asset_ptr_to_id(d):
-    classes = [surface_material.SurfaceMaterial, volume_material.VolumeMaterial, texture.Texture]
+    classes = [surface_material.SurfaceMaterial, volume_material.VolumeMaterial, texture.Texture,
+               levelset_3d.LevelSet3D]
     for key in d:
         for c in classes:
             if isinstance(d[key], c):

src/python/export_dynamics.cpp

@@ -12,12 +12,15 @@
 #include <taichi/dynamics/mpm2d/mpm.h>
 #include <taichi/dynamics/mpm2d/mpm_particle.h>
 #include <taichi/dynamics/simulation3d.h>
+#include <taichi/common/asset_manager.h>
 
 PYBIND11_MAKE_OPAQUE(std::vector<taichi::RenderParticle>);
 
 TC_NAMESPACE_BEGIN
 
 void export_dynamics(py::module &m) {
+    m.def("register_levelset3d", &AssetManager::insert_asset<LevelSet3D>);
+
     py::class_<Fluid::Particle>(m, "FluidParticle")
             .def(py::init<Vector2, Vector2>())
             .def_readwrite("position", &Fluid::Particle::position)

src/python/export_math.cpp

@@ -130,8 +130,8 @@ void export_math(py::module &m) {
     EXPORT_ARRAY_2D_OF(real, 1);
 
 #define EXPORT_ARRAY_3D_OF(T, C) \
-    py::class_<Array3D<real>>(m, "Array3D" #T)  \
-            .def(py::init<int, int, int>()) \
+    py::class_<Array3D<real>> PyArray3D##T(m, "Array3D" #T);  \
+        PyArray3D##T.def(py::init<int, int, int>()) \
             .def("get_width", &Array3D<T>::get_width) \
             .def("get_height", &Array3D<T>::get_height);
 
@@ -184,13 +184,15 @@ void export_math(py::module &m) {
             .def(py::init<>())
             .def("initialize", &DynamicLevelSet3D::initialize);
 
-    py::class_<LevelSet3D, Array3D<real>>(m, "LevelSet3D")
+    py::class_<LevelSet3D, std::shared_ptr<LevelSet3D>>(m, "LevelSet3D", PyArray3Dreal)
             .def(py::init<int, int, int, Vector3>())
             .def("get_width", &LevelSet3D::get_width)
             .def("get_height", &LevelSet3D::get_height)
             .def("get", &LevelSet3D::get_copy)
             .def("set", static_cast<void (LevelSet3D::*)(int, int, int, const real &)>(&LevelSet3D::set))
             .def("add_sphere", &LevelSet3D::add_sphere)
+            .def("add_plane", &LevelSet3D::add_plane)
+            .def("global_increase", &LevelSet3D::global_increase)
             .def("get_gradient", &LevelSet3D::get_gradient)
             .def("rasterize", &LevelSet3D::rasterize)
             .def("sample", static_cast<real(LevelSet3D::*)(real, real, real) const>(&LevelSet3D::sample))

src/python/export_visual.cpp

@@ -18,6 +18,7 @@
 #include <taichi/common/asset_manager.h>
 
 #include <taichi/geometry/factory.h>
+#include <taichi/math/levelset_3d.h>
 
 PYBIND11_MAKE_OPAQUE(std::vector<taichi::RenderParticle>);
 PYBIND11_MAKE_OPAQUE(std::vector<taichi::Triangle>);

src/texture/sources.cpp

@@ -11,6 +11,7 @@
 #include <taichi/visual/scene_geometry.h>
 #include <taichi/visualization/image_buffer.h>
 #include <taichi/math/array_3d.h>
+#include <taichi/math/levelset_3d.h>
 #include <taichi/common/asset_manager.h>
 
 TC_NAMESPACE_BEGIN
@@ -357,5 +358,24 @@ class MeshTexture : public Texture {
 
 TC_IMPLEMENTATION(Texture, MeshTexture, "mesh")
 
+class LevelSet3DTexture : public Texture {
+protected:
+    std::shared_ptr<LevelSet3D> levelset;
+    Vector2 bounds;
+public:
+    void initialize(const Config &config) override {
+        Texture::initialize(config);
+        levelset = AssetManager::get_asset<LevelSet3D>(config.get_int("levelset"));
+        bounds = config.get_vec2("bounds");
+    }
+
+    virtual Vector4 sample(const Vector3 &coord) const override {
+        real value = levelset->sample_relative_coord(coord);
+        return (bounds.x < value && value < bounds.y) ? Vector4(1) : Vector4(0);
+    }
+};
+
+TC_IMPLEMENTATION(Texture, LevelSet3DTexture, "levelset3d")
+
 TC_NAMESPACE_END