From e151d37f4bbcdc1568802abe27685842cc732dc2 Mon Sep 17 00:00:00 2001 From: fangy14 Date: Tue, 4 Apr 2017 11:00:08 +0800 Subject: [PATCH 1/4] Finish 2D dynamic level set --- include/taichi/math/dynamic_levelset_2d.cpp | 49 +++++++++++++++++++++ include/taichi/math/dynamic_levelset_2d.h | 35 +++++++++++++++ include/taichi/math/levelset_2d.cpp | 25 +++++++++-- include/taichi/math/levelset_2d.h | 3 ++ 4 files changed, 108 insertions(+), 4 deletions(-) create mode 100644 include/taichi/math/dynamic_levelset_2d.cpp create mode 100644 include/taichi/math/dynamic_levelset_2d.h diff --git a/include/taichi/math/dynamic_levelset_2d.cpp b/include/taichi/math/dynamic_levelset_2d.cpp new file mode 100644 index 0000000000000..cbe3335665dab --- /dev/null +++ b/include/taichi/math/dynamic_levelset_2d.cpp @@ -0,0 +1,49 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#include "dynamic_levelset_2d.h" + +TC_NAMESPACE_BEGIN + +void DynamicLevelSet2D::initialize(real _t0, real _t1, const LevelSet2D &_ls0, const LevelSet2D &_ls1) { + t0 = _t0; + t1 = _t1; + levelset0 = std::make_shared(_ls0); + levelset1 = std::make_shared(_ls1); +} + +Vector2 DynamicLevelSet2D::get_spatial_gradient(const Vector2 &pos, real t) { + 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); +} + +real DynamicLevelSet2D::get_temporal_gradient(const Vector2 &pos, real t) { + real l1 = levelset0->get(pos); + real l2 = levelset1->get(pos); + return (l2 - l1) / (t1 - t0); +} + +Array2D DynamicLevelSet2D::rasterize(int width, int height, real t) { + Array2D r0 = levelset0->rasterize(width, height); + Array2D r1 = levelset1->rasterize(width, height); + Array2D out(width, height); + for (auto &ind : Region2D(0, width, 0, height, Vector2(0.5f, 0.5f))) { + out[ind] = lerp((t - t0) / (t1 - t0), r0[ind], r1[ind]); + if (std::isnan(out[ind])) { + out[ind] = std::numeric_limits::infinity(); + } + } + return out; +} + + +TC_NAMESPACE_END diff --git a/include/taichi/math/dynamic_levelset_2d.h b/include/taichi/math/dynamic_levelset_2d.h new file mode 100644 index 0000000000000..690b0d2b685fa --- /dev/null +++ b/include/taichi/math/dynamic_levelset_2d.h @@ -0,0 +1,35 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#pragma once + + +#include +#include +#include + + +TC_NAMESPACE_BEGIN + +class DynamicLevelSet2D { + public: + real t0, t1; + std::shared_ptr levelset0, levelset1; + + void initialize(real _t0, real _t1, const LevelSet2D &_ls0, const LevelSet2D &_ls1); + + Vector2 get_spatial_gradient(const Vector2 &pos, real t); + + real get_temporal_gradient(const Vector2 &pos, real t); + + Array2D rasterize(int width, int height, real t); + }; + +TC_NAMESPACE_END + diff --git a/include/taichi/math/levelset_2d.cpp b/include/taichi/math/levelset_2d.cpp index 36b0c6a9f50a6..766279ddcaace 100644 --- a/include/taichi/math/levelset_2d.cpp +++ b/include/taichi/math/levelset_2d.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. @@ -16,7 +17,7 @@ void LevelSet2D::add_sphere(Vector2 center, real radius, bool inside_out) { for (auto &ind : get_region()) { Vector2 sample = ind.get_pos(); real dist = (inside_out ? -1 : 1) * (length(center - sample) - radius); - set(ind, std::min(get(ind), dist)); + set(ind, std::min(Array2D::get(ind), dist)); } } @@ -25,7 +26,7 @@ void LevelSet2D::add_polygon(std::vector polygon, bool inside_out) for (auto &ind : get_region()) { Vector2 p = ind.get_pos(); real dist = ((inside_polygon(p, polygon) ^ inside_out) ? -1 : 1) * (nearest_distance(p, polygon)); - set(ind, std::min(get(ind), dist)); + set(ind, std::min(Array2D::get(ind), dist)); } } @@ -39,8 +40,8 @@ Vector2 LevelSet2D::get_gradient(const Vector2 &pos) const const int y_i = clamp(int(y), 0, height - 2); const real x_r = x - x_i; const real y_r = y - y_i; - const real gx = lerp(y_r, get(x_i + 1, y_i) - get(x_i, y_i), get(x_i + 1, y_i + 1) - get(x_i, y_i + 1)); - const real gy = lerp(x_r, get(x_i, y_i + 1) - get(x_i, y_i), get(x_i + 1, y_i + 1) - get(x_i + 1, y_i)); + const real gx = lerp(y_r, Array2D::get(x_i + 1, y_i) - Array2D::get(x_i, y_i), Array2D::get(x_i + 1, y_i + 1) - Array2D::get(x_i, y_i + 1)); + const real gy = lerp(x_r, Array2D::get(x_i, y_i + 1) - Array2D::get(x_i, y_i), Array2D::get(x_i + 1, y_i + 1) - Array2D::get(x_i + 1, y_i)); return Vector2(gx, gy); } @@ -53,6 +54,22 @@ Vector2 LevelSet2D::get_normalized_gradient(const Vector2 &pos) const return normalize(gradient); } +real LevelSet2D::get(const Vector2 &pos) const +{ + assert_info(inside(pos), "LevelSet Gradient Query out of Bound! (" + std::to_string(pos.x) + ", " + std::to_string(pos.y) + ")"); + real x = pos.x, y = pos.y; + x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps); + y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps); + const int x_i = clamp(int(x), 0, width - 2); + const int y_i = clamp(int(y), 0, height - 2); + const real x_r = x - x_i; + const real y_r = y - y_i; + const real ly0 = lerp(x_r, Array2D::get(x_i, y_i), Array2D::get(x_i + 1, y_i)); + const real ly1 = lerp(x_r, Array2D::get(x_i, y_i + 1), Array2D::get(x_i + 1, y_i + 1)); + return lerp(y_r, ly0, ly1); +} + + Array2D LevelSet2D::rasterize(int width, int height) { for (auto &p : (*this)) { if (std::isnan(p)) { diff --git a/include/taichi/math/levelset_2d.h b/include/taichi/math/levelset_2d.h index 4485e8ba387f7..e9b6267031e26 100644 --- a/include/taichi/math/levelset_2d.h +++ b/include/taichi/math/levelset_2d.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. @@ -49,6 +50,8 @@ class LevelSet2D : public Array2D { Vector2 get_normalized_gradient(const Vector2 &pos) const; + real get(const Vector2 &pos) const; + static real fraction_outside(real phi_a, real phi_b) { return 1.0f - fraction_inside(phi_a, phi_b); } From f0f13cbae18b52d0feb4909c54c5842abf393b31 Mon Sep 17 00:00:00 2001 From: fangy14 Date: Tue, 4 Apr 2017 15:47:07 +0800 Subject: [PATCH 2/4] Finish 3D level set --- include/taichi/math/levelset_2d.cpp | 2 +- include/taichi/math/levelset_3d.cpp | 113 ++++++++++++++++++++++++++++ include/taichi/math/levelset_3d.h | 71 +++++++++++++++++ 3 files changed, 185 insertions(+), 1 deletion(-) create mode 100644 include/taichi/math/levelset_3d.cpp create mode 100644 include/taichi/math/levelset_3d.h diff --git a/include/taichi/math/levelset_2d.cpp b/include/taichi/math/levelset_2d.cpp index 766279ddcaace..400e5b5d0b438 100644 --- a/include/taichi/math/levelset_2d.cpp +++ b/include/taichi/math/levelset_2d.cpp @@ -56,7 +56,7 @@ Vector2 LevelSet2D::get_normalized_gradient(const Vector2 &pos) const real LevelSet2D::get(const Vector2 &pos) const { - assert_info(inside(pos), "LevelSet Gradient Query out of Bound! (" + std::to_string(pos.x) + ", " + std::to_string(pos.y) + ")"); + assert_info(inside(pos), "LevelSet Query out of Bound! (" + std::to_string(pos.x) + ", " + std::to_string(pos.y) + ")"); real x = pos.x, y = pos.y; x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps); y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps); diff --git a/include/taichi/math/levelset_3d.cpp b/include/taichi/math/levelset_3d.cpp new file mode 100644 index 0000000000000..7cd88801646a3 --- /dev/null +++ b/include/taichi/math/levelset_3d.cpp @@ -0,0 +1,113 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#include "levelset_3d.h" + +TC_NAMESPACE_BEGIN + + +void LevelSet3D::add_sphere(Vector3 center, real radius, bool inside_out) { + for (auto &ind : get_region()) { + Vector3 sample = ind.get_pos(); + real dist = (inside_out ? -1 : 1) * (length(center - sample) - radius); + set(ind, std::min(Array3D::get(ind), dist)); + } +} + +Vector3 LevelSet3D::get_gradient(const Vector3 &pos) const +{ + assert_info(inside(pos), "LevelSet Gradient Query out of Bound! (" + + std::to_string(pos.x) + ", " + + std::to_string(pos.y) + ", " + + std::to_string(pos.z) + ")"); + real x = pos.x, y = pos.y, z = pos.z; + x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps); + y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps); + z = clamp(z - storage_offset.z, 0.f, depth - 1.f - eps); + const int x_i = clamp(int(x), 0, width - 2); + const int y_i = clamp(int(y), 0, height - 2); + const int z_i = clamp(int(z), 0, depth - 2); + const real x_r = x - x_i; + const real y_r = y - y_i; + const real z_r = z - z_i; + const real gx = lerp(y_r, + lerp(z_r, Array3D::get(x_i + 1, y_i, z_i) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i + 1, y_i, z_i + 1) - Array3D::get(x_i, y_i, z_i + 1)), + lerp(z_r, Array3D::get(x_i + 1, y_i + 1, z_i) - Array3D::get(x_i, y_i + 1, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i + 1, z_i + 1))); + const real gy = lerp(z_r, + lerp(x_r, Array3D::get(x_i, y_i + 1, z_i) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i + 1, y_i + 1, z_i) - Array3D::get(x_i + 1, y_i, z_i)), + lerp(x_r, Array3D::get(x_i, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i, z_i + 1), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i + 1, y_i, z_i + 1))); + const real gz = lerp(x_r, + lerp(y_r, Array3D::get(x_i, y_i, z_i + 1) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i + 1, z_i)), + lerp(y_r, Array3D::get(x_i + 1, y_i, z_i + 1) - Array3D::get(x_i + 1, y_i, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i + 1, y_i + 1, z_i))); + return Vector3(gx, gy, gz); +} + +Vector3 LevelSet3D::get_normalized_gradient(const Vector3 &pos) const +{ + Vector3 gradient = get_gradient(pos); + if (length(gradient) < 1e-10f) + return Vector3(1, 0, 0); + else + return normalize(gradient); +} + +real LevelSet3D::get(const Vector3 &pos) const +{ + assert_info(inside(pos), "LevelSet Query out of Bound! (" + + std::to_string(pos.x) + ", " + + std::to_string(pos.y) + ", " + + std::to_string(pos.z) + ")"); + real x = pos.x, y = pos.y, z = pos.z; + x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps); + y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps); + z = clamp(z - storage_offset.z, 0.f, depth - 1.f - eps); + const int x_i = clamp(int(x), 0, width - 2); + const int y_i = clamp(int(y), 0, height - 2); + const int z_i = clamp(int(z), 0, depth - 2); + const real x_r = x - x_i; + const real y_r = y - y_i; + const real z_r = z - z_i; + return lerp(x_r, + lerp(y_r, + lerp(z_r, Array3D::get(x_i, y_i, z_i), Array3D::get(x_i, y_i, z_i + 1)), + lerp(z_r, Array3D::get(x_i, y_i + 1, z_i), Array3D::get(x_i, y_i + 1, z_i + 1))), + lerp(y_r, + lerp(z_r, Array3D::get(x_i + 1, y_i, z_i), Array3D::get(x_i + 1, y_i, z_i + 1)), + lerp(z_r, Array3D::get(x_i + 1, y_i + 1, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1)))); +} + + +Array3D LevelSet3D::rasterize(int width, int height, int depth) { + for (auto &p : (*this)) { + if (std::isnan(p)) { + printf("Warning: nan in levelset."); + } + } + Array3D out(width, height, depth); + Vector3 actual_size; + if (storage_offset == Vector3(0.0f, 0.0f, 0.0f)) { + actual_size = Vector3(this->width - 1, this->height - 1, this->depth - 1); + } + else { + actual_size = Vector3(this->width, this->height, this->depth); + } + + Vector3 scale_factor = actual_size / Vector3(width, height, depth); + + for (auto &ind : Region3D(0, width, 0, height, 0, depth, Vector3(0.5f, 0.5f, 0.5f))) { + Vector3 p = scale_factor * ind.get_pos(); + out[ind] = sample(p); + if (std::isnan(out[ind])) { + out[ind] = std::numeric_limits::infinity(); + } + } + return out; +} + +TC_NAMESPACE_END diff --git a/include/taichi/math/levelset_3d.h b/include/taichi/math/levelset_3d.h new file mode 100644 index 0000000000000..3b0ff02653c1f --- /dev/null +++ b/include/taichi/math/levelset_3d.h @@ -0,0 +1,71 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#pragma once + +#include +#include + + +TC_NAMESPACE_BEGIN + +class LevelSet3D : public Array3D { +public: + real friction = 1.0f; + + LevelSet3D() : LevelSet3D(0, 0, 0) {} + + LevelSet3D(int width, int height, int depth, Vector3 offset = Vector3(0.5f, 0.5f, 0.5f)) { + initialize(width, height, depth, offset); + } + + void initialize(int width, int height, int depth, Vector3 offset) { + Array3D::initialize(width, height, depth, std::numeric_limits::infinity(), offset); + } + + void initialize(const Vector3i &res, Vector3 offset) { + Array3D::initialize(res[0], res[1], res[2], std::numeric_limits::infinity(), offset); + } + + void initialize(int width, int height, int depth, Vector3 offset, real value) { + Array3D::initialize(width, height, depth, value, offset); + } + + void initialize(const Vector3i &res, Vector3 offset, real value) { + Array3D::initialize(res[0], res[1], res[2], value, offset); + } + + void add_sphere(Vector3 center, real radius, bool inside_out = false); + + Vector3 get_gradient(const Vector3 &pos) const; // Note this is not normalized! + + Vector3 get_normalized_gradient(const Vector3 &pos) const; + + real get(const Vector3 &pos) const; + + static real fraction_outside(real phi_a, real phi_b) { + return 1.0f - fraction_inside(phi_a, phi_b); + } + + static real fraction_inside(real phi_a, real phi_b) { + if (phi_a < 0 && phi_b < 0) + return 1; + if (phi_a < 0 && phi_b >= 0) + return phi_a / (phi_a - phi_b); + if (phi_a >= 0 && phi_b < 0) + return phi_b / (phi_b - phi_a); + else + return 0; + } + + Array3D rasterize(int width, int height, int depth); +}; + +TC_NAMESPACE_END + From 4b05f705af1f2679f6b724fc2dc38765e5f841a4 Mon Sep 17 00:00:00 2001 From: fangy14 Date: Tue, 4 Apr 2017 16:06:04 +0800 Subject: [PATCH 3/4] Finish 3D dynamic level set --- include/taichi/math/dynamic_levelset_3d.cpp | 50 +++++++++++++++++++++ include/taichi/math/dynamic_levelset_3d.h | 35 +++++++++++++++ 2 files changed, 85 insertions(+) create mode 100644 include/taichi/math/dynamic_levelset_3d.cpp create mode 100644 include/taichi/math/dynamic_levelset_3d.h diff --git a/include/taichi/math/dynamic_levelset_3d.cpp b/include/taichi/math/dynamic_levelset_3d.cpp new file mode 100644 index 0000000000000..cd346b313a0dc --- /dev/null +++ b/include/taichi/math/dynamic_levelset_3d.cpp @@ -0,0 +1,50 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#include "dynamic_levelset_3d.h" + +TC_NAMESPACE_BEGIN + +void DynamicLevelSet3D::initialize(real _t0, real _t1, const LevelSet3D &_ls0, const LevelSet3D &_ls1) { + t0 = _t0; + t1 = _t1; + levelset0 = std::make_shared(_ls0); + levelset1 = std::make_shared(_ls1); +} + +Vector3 DynamicLevelSet3D::get_spatial_gradient(const Vector3 &pos, real t) { + 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); +} + +real DynamicLevelSet3D::get_temporal_gradient(const Vector3 &pos, real t) { + real l1 = levelset0->get(pos); + real l2 = levelset1->get(pos); + return (l2 - l1) / (t1 - t0); +} + +Array3D DynamicLevelSet3D::rasterize(int width, int height, int depth, real t) { + Array3D r0 = levelset0->rasterize(width, height, depth); + Array3D r1 = levelset1->rasterize(width, height, depth); + Array3D out(width, height, width); + 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])) { + out[ind] = std::numeric_limits::infinity(); + } + } + return out; +} + + +TC_NAMESPACE_END diff --git a/include/taichi/math/dynamic_levelset_3d.h b/include/taichi/math/dynamic_levelset_3d.h new file mode 100644 index 0000000000000..73c2a8acace9c --- /dev/null +++ b/include/taichi/math/dynamic_levelset_3d.h @@ -0,0 +1,35 @@ +/******************************************************************************* + Taichi - Physically based Computer Graphics Library + + Copyright (c) 2017 Yu Fang + + All rights reserved. Use of this source code is governed by + the MIT license as written in the LICENSE file. +*******************************************************************************/ + +#pragma once + + +#include +#include +#include + + +TC_NAMESPACE_BEGIN + +class DynamicLevelSet3D { +public: + real t0, t1; + std::shared_ptr levelset0, levelset1; + + void initialize(real _t0, real _t1, const LevelSet3D &_ls0, const LevelSet3D &_ls1); + + Vector3 get_spatial_gradient(const Vector3 &pos, real t); + + real get_temporal_gradient(const Vector3 &pos, real t); + + Array3D rasterize(int width, int height, int depth, real t); +}; + +TC_NAMESPACE_END + From 20f19dfd30c19e7bf558fc7d88aa8b63ad92043d Mon Sep 17 00:00:00 2001 From: fangy14 Date: Thu, 6 Apr 2017 14:14:10 +0800 Subject: [PATCH 4/4] Replace static levelset with dynamic levelset for 2D mpm simulator --- include/taichi/dynamics/mpm2d/mpm.cpp | 7 +++-- include/taichi/dynamics/mpm2d/mpm.h | 6 ++-- include/taichi/dynamics/mpm2d/mpm_grid.cpp | 13 +++++---- include/taichi/dynamics/mpm2d/mpm_grid.h | 4 ++- include/taichi/dynamics/mpm2d/mpm_particle.h | 8 ++++-- include/taichi/math/dynamic_levelset_2d.cpp | 20 +++++++++---- include/taichi/math/dynamic_levelset_2d.h | 8 ++++-- include/taichi/math/dynamic_levelset_3d.cpp | 18 +++++++++--- include/taichi/math/dynamic_levelset_3d.h | 6 ++-- .../examples/simulation/2d/mpm_sand_stir.py | 28 +++++++++++++++++++ python/taichi/two_d/mpm_simulator.py | 6 ++-- python/taichi/two_d/simulator.py | 21 +++++++++++--- src/python/export_math.cpp | 5 ++++ 13 files changed, 116 insertions(+), 34 deletions(-) create mode 100644 python/examples/simulation/2d/mpm_sand_stir.py 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);