diff --git a/include/taichi/dynamics/mpm2d/mpm.cpp b/include/taichi/dynamics/mpm2d/mpm.cpp index c5f4b28f79363..f0532dd75f5f8 100644 --- a/include/taichi/dynamics/mpm2d/mpm.cpp +++ b/include/taichi/dynamics/mpm2d/mpm.cpp @@ -2,6 +2,7 @@ Taichi - Physically based Computer Graphics Library Copyright (c) 2016 Yuanming Hu + 2017 Yu Fang 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); } } diff --git a/include/taichi/dynamics/mpm2d/mpm.h b/include/taichi/dynamics/mpm2d/mpm.h index ffb73ebf244d7..dd87f244b740f 100644 --- a/include/taichi/dynamics/mpm2d/mpm.h +++ b/include/taichi/dynamics/mpm2d/mpm.h @@ -2,6 +2,7 @@ Taichi - Physically based Computer Graphics Library Copyright (c) 2016 Yuanming Hu + 2017 Yu Fang 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 #include "mpm_grid.h" #include +#include #include #include @@ -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; } diff --git a/include/taichi/dynamics/mpm2d/mpm_grid.cpp b/include/taichi/dynamics/mpm2d/mpm_grid.cpp index 0a9d61a992e72..56ef5e23af6c0 100644 --- a/include/taichi/dynamics/mpm2d/mpm_grid.cpp +++ b/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; } } diff --git a/include/taichi/dynamics/mpm2d/mpm_grid.h b/include/taichi/dynamics/mpm2d/mpm_grid.h index c970ee5f74551..4a2b151129b82 100644 --- a/include/taichi/dynamics/mpm2d/mpm_grid.h +++ b/include/taichi/dynamics/mpm2d/mpm_grid.h @@ -2,6 +2,7 @@ Taichi - Physically based Computer Graphics Library Copyright (c) 2016 Yuanming Hu + 2017 Yu Fang 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 #include #include +#include #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++) { diff --git a/include/taichi/dynamics/mpm2d/mpm_particle.h b/include/taichi/dynamics/mpm2d/mpm_particle.h index 661585c79a631..cc92bd42b2250 100644 --- a/include/taichi/dynamics/mpm2d/mpm_particle.h +++ b/include/taichi/dynamics/mpm2d/mpm_particle.h @@ -2,6 +2,7 @@ Taichi - Physically based Computer Graphics Library Copyright (c) 2016 Yuanming Hu + 2017 Yu Fang 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 #include +#include 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; } } diff --git a/include/taichi/math/dynamic_levelset_2d.cpp b/include/taichi/math/dynamic_levelset_2d.cpp index cbe3335665dab..d183082561866 100644 --- a/include/taichi/math/dynamic_levelset_2d.cpp +++ b/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(_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 DynamicLevelSet2D::rasterize(int width, int height, real t) { diff --git a/include/taichi/math/dynamic_levelset_2d.h b/include/taichi/math/dynamic_levelset_2d.h index 690b0d2b685fa..b1bfa47dfd8d9 100644 --- a/include/taichi/math/dynamic_levelset_2d.h +++ b/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 rasterize(int width, int height, real t); - }; +}; TC_NAMESPACE_END diff --git a/include/taichi/math/dynamic_levelset_3d.cpp b/include/taichi/math/dynamic_levelset_3d.cpp index cd346b313a0dc..96bc9c4463f5c 100644 --- a/include/taichi/math/dynamic_levelset_3d.cpp +++ b/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(_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 DynamicLevelSet3D::rasterize(int width, int height, int depth, real t) { diff --git a/include/taichi/math/dynamic_levelset_3d.h b/include/taichi/math/dynamic_levelset_3d.h index 73c2a8acace9c..c1f9b043c6df8 100644 --- a/include/taichi/math/dynamic_levelset_3d.h +++ b/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 rasterize(int width, int height, int depth, real t); }; diff --git a/python/examples/simulation/2d/mpm_sand_stir.py b/python/examples/simulation/2d/mpm_sand_stir.py new file mode 100644 index 0000000000000..a7ee2e313a46c --- /dev/null +++ b/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) diff --git a/python/taichi/two_d/mpm_simulator.py b/python/taichi/two_d/mpm_simulator.py index a59424bb9838a..659616d02dc78 100644 --- a/python/taichi/two_d/mpm_simulator.py +++ b/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'])) diff --git a/python/taichi/two_d/simulator.py b/python/taichi/two_d/simulator.py index 6542b9849a603..f4cb4d105cb30 100644 --- a/python/taichi/two_d/simulator.py +++ b/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): diff --git a/src/python/export_math.cpp b/src/python/export_math.cpp index b3cd367d3d39a..c02991ddf9eb3 100644 --- a/src/python/export_math.cpp +++ b/src/python/export_math.cpp @@ -12,6 +12,7 @@ #include #include #include +#include PYBIND11_MAKE_OPAQUE(std::vector); PYBIND11_MAKE_OPAQUE(std::vector); @@ -111,6 +112,10 @@ void export_math(py::module &m) { .def("to_ndarray", &array2d_to_ndarray) .def_readwrite("friction", &LevelSet2D::friction); + py::class_(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);