Add dummy levelset

include/taichi/dynamics/mpm2d/mpm_grid.cpp

@@ -14,34 +14,32 @@ 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) {
-        for (auto &ind : scheduler->get_active_grid_points()) {
-            Vector2 pos = Vector2(ind[0] + 0.5f, ind[1] + 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.levelset0->friction;
-                if (mu < 0) { // sticky
-                    v = Vector2(0.0f);
-                }
-                else {
-                    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;
-                }
-            }
-            else if (phi <= 0) {
+    for (auto &ind : scheduler->get_active_grid_points()) {
+        Vector2 pos = Vector2(ind[0] + 0.5f, ind[1] + 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.levelset0->friction;
+            if (mu < 0) { // sticky
                 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;
+            else {
+                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;
+            }
+        }
+        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/scheduler.h

@@ -212,11 +212,12 @@ class MPMScheduler {
             for (int i = 0; i < 4; i++) {
                 block_absolute_vel = std::max(block_absolute_vel, std::abs(min_max_vel_expanded[ind][i]));
             }
-            real distance2boundary = std::max(
-                    levelset->sample(Vector2(ind.get_pos() * real(grid_block_size)), t) - real(grid_block_size) * 0.75f,
-                    0.5f);
-            int64 boundary_limit = int64(cfl * distance2boundary / block_absolute_vel / base_delta_t);
-            cfl_limit = std::min(cfl_limit, boundary_limit);
+			real last_distance = levelset->sample(Vector2(ind.get_pos() * real(grid_block_size)), t);
+			if (last_distance < LevelSet2D::INF) {
+				real distance2boundary = std::max(last_distance - real(grid_block_size) * 0.75f, 0.5f);
+				int64 boundary_limit = int64(cfl * distance2boundary / block_absolute_vel / base_delta_t);
+				cfl_limit = std::min(cfl_limit, boundary_limit);
+			}
             max_dt_int_cfl[ind] = get_largest_pot(cfl_limit);
         }
 	}

include/taichi/math/levelset_2d.cpp

@@ -12,6 +12,7 @@
 
 TC_NAMESPACE_BEGIN
 
+const real LevelSet2D::INF = 1e7;
 
 void LevelSet2D::add_sphere(Vector2 center, real radius, bool inside_out) {
     for (auto &ind : get_region()) {

include/taichi/math/levelset_2d.h

@@ -18,6 +18,8 @@ TC_NAMESPACE_BEGIN
 
 class LevelSet2D : public Array2D<real> {
 public:
+    static const real INF;
+
     real friction = 1.0f;
 
     LevelSet2D() : LevelSet2D(0, 0) {}
@@ -27,11 +29,11 @@ class LevelSet2D : public Array2D<real> {
     }
 
     void initialize(int width, int height, Vector2 offset) {
-        Array2D<real>::initialize(width, height, std::numeric_limits<real>::infinity(), offset);
+        Array2D<real>::initialize(width, height, INF, offset);
     }
 
     void initialize(const Vector2i &res, Vector2 offset) {
-        Array2D<real>::initialize(res[0], res[1], std::numeric_limits<real>::infinity(), offset);
+        Array2D<real>::initialize(res[0], res[1], INF, offset);
     }
 
     void initialize(int width, int height, Vector2 offset, real value) {

include/taichi/math/levelset_3d.cpp

@@ -11,6 +11,7 @@
 
 TC_NAMESPACE_BEGIN
 
+const real LevelSet3D::INF = 1e7;
 
 void LevelSet3D::add_sphere(Vector3 center, real radius, bool inside_out) {
     for (auto &ind : get_region()) {

include/taichi/math/levelset_3d.h

@@ -17,6 +17,8 @@ TC_NAMESPACE_BEGIN
 
 class LevelSet3D : public Array3D<real> {
 public:
+    static const real INF;
+
     real friction = 1.0f;
 
     LevelSet3D() : LevelSet3D(0, 0, 0) {}
@@ -26,11 +28,11 @@ class LevelSet3D : public Array3D<real> {
     }
 
     void initialize(int width, int height, int depth, Vector3 offset) {
-        Array3D<real>::initialize(width, height, depth, std::numeric_limits<real>::infinity(), offset);
+        Array3D<real>::initialize(width, height, depth, INF, offset);
     }
 
     void initialize(const Vector3i &res, Vector3 offset) {
-        Array3D<real>::initialize(res[0], res[1], res[2], std::numeric_limits<real>::infinity(), offset);
+        Array3D<real>::initialize(res[0], res[1], res[2], INF, offset);
     }
 
     void initialize(int width, int height, int depth, Vector3 offset, real value) {

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

@@ -8,11 +8,11 @@
 from taichi.visual.texture import Texture
 import taichi as tc
 
-gi_render = False
+gi_render = True
 step_number = 400
 # step_number = 1
 # total_frames = 1
-grid_downsample = 8
+grid_downsample = 3
 output_downsample = 1
 render_epoch = 20
 
@@ -53,7 +53,7 @@ def create_sand_scene(frame, d, t):
 
         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))
+                         translate=((t+0.05) * 0.5 - 0.35, -0.60, 0), scale=0.1, rotation=(0, 0, 0))
         scene.add_mesh(sphere)
 
         # Change this line to your particle output path pls.
@@ -75,10 +75,10 @@ def render_sand_frame(frame, d, t):
 if __name__ == '__main__':
     downsample = grid_downsample
     resolution = (255 / downsample, 255 / downsample, 255 / downsample)
-    tex = Texture('ring', outer=0.15) * 4
+    tex = Texture('ring', outer=0.15) * 8
     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,
+    mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), initial_velocity=(0, 0, 0), delta_t=0.001, num_threads=8,
                density_tex=tex.id)
 
 

python/taichi/dynamics/levelset_3d.py

@@ -5,7 +5,7 @@
 class LevelSet3D:
     def __init__(self, res, offset=None):
         if offset is None:
-            offset = Vector(0, 0, 0)
+            offset = Vector(0.5, 0.5, 0.5)
         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)

python/taichi/dynamics/mpm.py

@@ -28,7 +28,11 @@ def __init__(self, **kwargs):
                                                   light_direction=(1, 1, 0))
         self.resolution = kwargs['resolution']
         self.frame = 0
-        self.levelset_generator = None
+
+        def dummy_levelset_generator(_):
+            return self.create_levelset()
+
+        self.levelset_generator = dummy_levelset_generator
 
     def update_levelset(self, t0, t1):
         levelset = tc.core.DynamicLevelSet3D()

python/taichi/two_d/mpm_simulator.py

@@ -18,6 +18,11 @@ def __init__(self, **kwargs):
         self.delta_x = kwargs['delta_x']
         self.sample_rate = kwargs.get('sample_rate', 2)
 
+        def dummy_levelset_generator(_):
+            return self.create_levelset()
+
+        self.levelset_generator = dummy_levelset_generator
+
     @staticmethod
     def create_particle(particle_type):
         if particle_type == 'ep':

src/simulation3d/mpm/mpm3.cpp

@@ -312,42 +312,38 @@ void MPM3D::apply_deformation_force(float delta_t) {
 }
 
 void MPM3D::grid_apply_boundary_conditions(const DynamicLevelSet3D &levelset, real t) {
-    if (levelset.levelset0) {
-        for (auto &ind : grid_velocity.get_region()) {
-            Vector3 pos = Vector3(ind.get_pos());
-            Vector3 v = grid_velocity[ind] + grid_force_or_acc[ind] * delta_t -
-                        levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
-            Vector3 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.levelset0->friction;
-                if (mu < 0) { // sticky
-                    v = Vector3(0.0f);
-                } else {
-                    Vector3 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;
-                }
-            } else if (phi <= 0) {
+    for (auto &ind : grid_velocity.get_region()) {
+        Vector3 pos = Vector3(ind.get_pos());
+        Vector3 v = grid_velocity[ind] + grid_force_or_acc[ind] * delta_t -
+                    levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
+        Vector3 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.levelset0->friction;
+            if (mu < 0) { // sticky
                 v = Vector3(0.0f);
+            } else {
+                Vector3 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;
             }
-            v += levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
-            grid_force_or_acc[ind] = (v - grid_velocity[ind]) / delta_t;
+        } else if (phi <= 0) {
+            v = Vector3(0.0f);
         }
+        v += levelset.get_temporal_derivative(pos, t) * levelset.get_spatial_gradient(pos, t);
+        grid_force_or_acc[ind] = (v - grid_velocity[ind]) / delta_t;
     }
 }
 
 void MPM3D::particle_collision_resolution(real t) {
-    if (levelset.levelset0) {
-        parallel_for_each_particle([&](Particle &p) {
-            p.resolve_collision(levelset, t);
-        });
-    }
+    parallel_for_each_particle([&](Particle &p) {
+        p.resolve_collision(levelset, t);
+    });
 }
 
 void MPM3D::substep(float delta_t) {