Skip to content

Commit

Permalink
Merge pull request #37 from squarefk/ampm
Browse files Browse the repository at this point in the history
2D/3D level set/dynamic level set for MPM
  • Loading branch information
yuanming-hu committed Apr 6, 2017
2 parents 1d8c85b + 20f19df commit d323a57
Show file tree
Hide file tree
Showing 17 changed files with 479 additions and 24 deletions.
7 changes: 4 additions & 3 deletions include/taichi/dynamics/mpm2d/mpm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
Taichi - Physically based Computer Graphics Library
Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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);
}
}

Expand Down
6 changes: 4 additions & 2 deletions include/taichi/dynamics/mpm2d/mpm.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
Taichi - Physically based Computer Graphics Library
Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
Expand All @@ -13,6 +14,7 @@
#include <vector>
#include "mpm_grid.h"
#include <taichi/math/levelset_2d.h>
#include <taichi/math/dynamic_levelset_2d.h>
#include <taichi/visual/texture.h>
#include <taichi/visualization/image_buffer.h>

Expand Down Expand Up @@ -42,7 +44,7 @@ class MPM {
Vector2 gravity;
bool apic;

LevelSet2D levelset;
DynamicLevelSet2D levelset;
LevelSet2D material_levelset;

real last_sort;
Expand Down Expand Up @@ -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;
}

Expand Down
13 changes: 8 additions & 5 deletions include/taichi/dynamics/mpm2d/mpm_grid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand All @@ -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;
}
}
Expand Down
4 changes: 3 additions & 1 deletion include/taichi/dynamics/mpm2d/mpm_grid.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
Taichi - Physically based Computer Graphics Library
Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
Expand All @@ -16,6 +17,7 @@
#include <taichi/math/array_2d.h>
#include <taichi/math/array_1d.h>
#include <taichi/math/levelset_2d.h>
#include <taichi/math/dynamic_levelset_2d.h>
#include "mpm_particle.h"

TC_NAMESPACE_BEGIN
Expand Down Expand Up @@ -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++) {
Expand Down
8 changes: 5 additions & 3 deletions include/taichi/dynamics/mpm2d/mpm_particle.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
Taichi - Physically based Computer Graphics Library
Copyright (c) 2016 Yuanming Hu <yuanmhu@gmail.com>
2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
Expand All @@ -14,6 +15,7 @@
#include "mpm_utils.h"
#include <taichi/common/util.h>
#include <taichi/math/levelset_2d.h>
#include <taichi/math/dynamic_levelset_2d.h>

TC_NAMESPACE_BEGIN

Expand Down Expand Up @@ -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;
}
}

Expand Down
59 changes: 59 additions & 0 deletions include/taichi/math/dynamic_levelset_2d.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
/*******************************************************************************
Taichi - Physically based Computer Graphics Library
Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
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<LevelSet2D>(_ls0);
levelset1 = std::make_shared<LevelSet2D>(_ls1);
}

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);
Vector2 gradient = Vector2(gx, gy);
if (length(gradient) < 1e-10f)
return Vector2(1, 0);
else
return normalize(gradient);
}

real DynamicLevelSet2D::get_temporal_derivative(const Vector2 &pos, real t) const {
real l1 = levelset0->get(pos);
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<real> DynamicLevelSet2D::rasterize(int width, int height, real t) {
Array2D<real> r0 = levelset0->rasterize(width, height);
Array2D<real> r1 = levelset1->rasterize(width, height);
Array2D<real> 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<real>::infinity();
}
}
return out;
}


TC_NAMESPACE_END
37 changes: 37 additions & 0 deletions include/taichi/math/dynamic_levelset_2d.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
/*******************************************************************************
Taichi - Physically based Computer Graphics Library
Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
*******************************************************************************/

#pragma once


#include <limits>
#include <memory>
#include <taichi/math/levelset_2d.h>


TC_NAMESPACE_BEGIN

class DynamicLevelSet2D {
public:
real t0, t1;
std::shared_ptr<LevelSet2D> levelset0, levelset1;

void initialize(real _t0, real _t1, const LevelSet2D &_ls0, const LevelSet2D &_ls1);

Vector2 get_spatial_gradient(const Vector2 &pos, real t) const;

real get_temporal_derivative(const Vector2 &pos, real t) const;

real sample(const Vector2 &pos, real t) const;

Array2D<real> rasterize(int width, int height, real t);
};

TC_NAMESPACE_END

60 changes: 60 additions & 0 deletions include/taichi/math/dynamic_levelset_3d.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
/*******************************************************************************
Taichi - Physically based Computer Graphics Library
Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
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<LevelSet3D>(_ls0);
levelset1 = std::make_shared<LevelSet3D>(_ls1);
}

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);
Vector3 gradient = Vector3(gx, gy, gz);
if (length(gradient) < 1e-10f)
return Vector3(1, 0, 0);
else
return normalize(gradient);
}

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 lerp((t - t0) / (t1 - t0), l1, l2);
}

Array3D<real> DynamicLevelSet3D::rasterize(int width, int height, int depth, real t) {
Array3D<real> r0 = levelset0->rasterize(width, height, depth);
Array3D<real> r1 = levelset1->rasterize(width, height, depth);
Array3D<real> 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<real>::infinity();
}
}
return out;
}


TC_NAMESPACE_END
37 changes: 37 additions & 0 deletions include/taichi/math/dynamic_levelset_3d.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
/*******************************************************************************
Taichi - Physically based Computer Graphics Library
Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
All rights reserved. Use of this source code is governed by
the MIT license as written in the LICENSE file.
*******************************************************************************/

#pragma once


#include <limits>
#include <memory>
#include <taichi/math/levelset_3d.h>


TC_NAMESPACE_BEGIN

class DynamicLevelSet3D {
public:
real t0, t1;
std::shared_ptr<LevelSet3D> levelset0, levelset1;

void initialize(real _t0, real _t1, const LevelSet3D &_ls0, const LevelSet3D &_ls1);

Vector3 get_spatial_gradient(const Vector3 &pos, real t) const;

real get_temporal_derivative(const Vector3 &pos, real t) const;

real sample(const Vector3 &pos, real t) const;

Array3D<real> rasterize(int width, int height, int depth, real t);
};

TC_NAMESPACE_END

Loading

0 comments on commit d323a57

Please sign in to comment.