Merge pull request #42 from squarefk/ampm

Ampm

.travis.yml

@@ -25,8 +25,10 @@ before_script:
   - export LD_LIBRARY_PATH=/usr/lib64:$LD_LIBRARY_PATH
   - mkdir build
   - cd build
-  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then sudo sh ../build_scripts/linux.sh ; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then sudo sh ../build_scripts/osx.sh ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then travis_wait 20 sudo sh ../scripts/linux_sudo.sh ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then sudo sh ../scripts/osx_sudo.sh ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then sh ../scripts/linux_normal.sh ; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then sh ../scripts/osx_normal.sh ; fi
 
 script:
   - python -c 'import taichi as tc'

CMakeLists.txt

@@ -96,7 +96,7 @@ file(GLOB TAICHI_SOURCE
         "include/taichi/*/*/*/*.h" "include/taichi/*/*/*.h" "include/taichi/*/*.h" "include/taichi/*.h")
 
 set(CORE_LIBRARY_NAME taichi_core)
-add_library(${CORE_LIBRARY_NAME} SHARED ${TAICHI_SOURCE} include/taichi/dynamics/mpm2d/scheduler.cpp)
+add_library(${CORE_LIBRARY_NAME} SHARED ${TAICHI_SOURCE})
 
 if (NOT WIN32)
     target_link_libraries(${CORE_LIBRARY_NAME} pthread)

include/taichi/dynamics/mpm2d/mpm.h

@@ -12,7 +12,7 @@
 
 #include <memory>
 #include <vector>
-#include "scheduler.h"
+#include "mpm_scheduler.h"
 #include "mpm_grid.h"
 #include <taichi/math/levelset_2d.h>
 #include <taichi/math/dynamic_levelset_2d.h>

include/taichi/dynamics/mpm2d/mpm_grid.h

@@ -11,15 +11,14 @@
 #pragma once
 
 #include "mpm_utils.h"
-#include <stb_image.h>
 #include <algorithm>
 #include <atomic>
 #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"
-#include "scheduler.h"
+#include "mpm_scheduler.h"
 
 TC_NAMESPACE_BEGIN
 

include/taichi/dynamics/mpm2d/mpm_particle.h

@@ -14,6 +14,7 @@
 #include <taichi/math/array_2d.h>
 #include "mpm_utils.h"
 #include <taichi/common/util.h>
+#include <taichi/math/qr_svd.h>
 #include <taichi/math/levelset_2d.h>
 #include <taichi/math/dynamic_levelset_2d.h>
 

include/taichi/dynamics/mpm2d/scheduler.cpp → include/taichi/dynamics/mpm2d/mpm_scheduler.cpp

@@ -7,11 +7,11 @@
     the MIT license as written in the LICENSE file.
 *******************************************************************************/
 
-#include "scheduler.h"
+#include "mpm_scheduler.h"
 
 TC_NAMESPACE_BEGIN
 
-template <typename T> using Array = Array2D<T>;
+template<typename T> using Array = Array2D<T>;
 
 void MPMScheduler::expand(bool expand_vel, bool expand_state) {
     Array<int> new_states;

include/taichi/dynamics/mpm2d/scheduler.h → include/taichi/dynamics/mpm2d/mpm_scheduler.h

@@ -20,7 +20,7 @@ class MPMScheduler {
 public:
     typedef MPMParticle Particle;
 
-    template <typename T> using Array = Array2D<T>;
+    template<typename T> using Array = Array2D<T>;
 
     Array<int64> max_dt_int_strength;
     Array<int64> max_dt_int_cfl;

include/taichi/dynamics/mpm2d/mpm_utils.h

@@ -65,15 +65,5 @@ inline Vector2 dw(const Vector2 &a) {
     return Vector2(dw(a.x) * w(a.y), w(a.x) * dw(a.y));
 }
 
-inline real det(const Matrix2 &m) {
-    return m[0][0] * m[1][1] - m[0][1] * m[1][0];
-}
-
-void polar_decomp(const Matrix2 &A, Matrix2 &r, Matrix2 &s);
-
-void svd(Matrix2 A, Matrix2 &u, Matrix2 &sig, Matrix2 &v);
-void svd(Matrix3 A, Matrix3 &u, Matrix3 &sig, Matrix3 &v);
-
-
 TC_NAMESPACE_END
 

include/taichi/dynamics/simulation3d.h

@@ -33,6 +33,10 @@ class Simulation3D : public Unit {
         num_threads = config.get_int("num_threads");
     }
 
+    virtual void add_particles(const Config &config) {
+        error("no impl");
+    }
+
     virtual void step(real t) {
         error("no impl");
     }

include/taichi/math/math_util.h

@@ -280,6 +280,14 @@ inline double pow(const double &a, const double &b) {
     return ::pow(a, b);
 }
 
+inline real det(const Matrix2 &m) {
+    return glm::determinant(m);
+}
+
+inline real det(const Matrix3 &m) {
+    return glm::determinant(m);
+}
+
 //#define rand frand
 
 

include/taichi/math/qr_svd.cpp

@@ -14,12 +14,6 @@
 TC_NAMESPACE_BEGIN
 
 void imp_svd(Matrix3 m, Matrix3 &u, Matrix3 &s, Matrix3 &v) {
-    /*
-    Eigen::Matrix<T, 3, 3> *M;
-    Eigen::Matrix<T, 3, 1> *S;
-    Eigen::Matrix<T, 3, 3> *U;
-    Eigen::Matrix<T, 3, 3> *V;
-    */
     JIXIE::singularValueDecomposition(
             *(Eigen::Matrix<float, 3, 3> *)&m,
             *(Eigen::Matrix<float, 3, 3> *)&u,
@@ -30,40 +24,17 @@ void imp_svd(Matrix3 m, Matrix3 &u, Matrix3 &s, Matrix3 &v) {
     memset(&s[0][0] + 1, 0, sizeof(float) * 8);
     s[1][1] = s_tmp[0];
     s[2][2] = s_tmp[1];
-    /*
-    printf("glm\n");
-    for (int i = 0; i < 3; i++) {
-        for (int j = 0; j < 3; j++) {
-            printf(" %.4f", m[j][i]);
-        }
-        printf("\n");
-    }
-    printf("\n");
-    printf("eigen\n");
-    for (int i = 0; i < 3; i++) {
-        for (int j = 0; j < 3; j++) {
-            printf(" %.4f", (*((Eigen::Matrix<float, 3, 3> *) &m))(i, j));
-        }
-        printf("\n");
-    }
-     */
 }
 
 // m can not be const here, otherwise JIXIE::singularValueDecomposition will cause a error due to const_cast
 void imp_svd(Matrix2 m, Matrix2 &u, Matrix2 &s, Matrix2 &v) {
-    /*
-    Eigen::Matrix<T, 3, 3> *M;
-    Eigen::Matrix<T, 3, 1> *S;
-    Eigen::Matrix<T, 3, 3> *U;
-    Eigen::Matrix<T, 3, 3> *V;
-    */
     JIXIE::singularValueDecomposition(
             *(Eigen::Matrix<float, 2, 2> *)&m,
             *(Eigen::Matrix<float, 2, 2> *)&u,
             *(Eigen::Matrix<float, 2, 1> *)&s,
             *(Eigen::Matrix<float, 2, 2> *)&v
     );
-    float s_tmp {s[0][1]};
+    float s_tmp{s[0][1]};
     memset(&s[0][0] + 1, 0, sizeof(float) * 3);
     if (s_tmp > 0) {
         s[1][1] = s_tmp;
@@ -74,4 +45,55 @@ void imp_svd(Matrix2 m, Matrix2 &u, Matrix2 &s, Matrix2 &v) {
     }
 }
 
+void svd(Matrix2 m, Matrix2 &u, Matrix2 &sig, Matrix2 &v) {
+    // TODO: what's going on ???
+    if (frobenius_norm2(m - Matrix2(m[0][0])) < 1e-7f) {
+        u = m;
+        sig = v = Matrix2(1);
+    } else {
+        imp_svd(m, u, sig, v);
+    }
+}
+
+void svd(Matrix3 A, Matrix3 &u, Matrix3 &sig, Matrix3 &v) {
+    if (frobenius_norm2(A - Matrix3(A[0][0])) < 1e-5f) {
+        u = A;
+        sig = v = Matrix3(1);
+    } else {
+        imp_svd(A, u, sig, v);
+    }
+}
+
+void polar_decomp(Matrix2 A, Matrix2 &r, Matrix2 &s) {
+    Matrix2 u, sig, v;
+    svd(A, u, sig, v);
+    r = u * glm::transpose(v);
+    s = v * sig * glm::transpose(v);
+}
+
+void polar_decomp(Matrix3 A, Matrix3 &r, Matrix3 &s) {
+    Matrix3 u, sig, v;
+    svd(A, u, sig, v);
+    r = u * glm::transpose(v);
+    s = v * sig * glm::transpose(v);
+    if (!is_normal(r)) {
+        Matrix3 m = A;
+        svd(m, u, sig, v);
+        P(A);
+        P(m);
+        P(u);
+        P(sig);
+        P(v);
+        P(r);
+        P(s);
+        P(glm::transpose(v));
+        P(u * glm::transpose(v));
+        r = u * glm::transpose(v);
+        P(r);
+        printf("Matrix3 m(%.30f,%.30f,%.30f,%.30f,%.30f,%.30f,%.30f,%.30f,%.30f);\n", m[0][0], m[1][0], m[2][0],
+               m[0][1],
+               m[1][1], m[2][1], m[0][2], m[1][2], m[2][2]);
+    }
+}
+
 TC_NAMESPACE_END

include/taichi/math/qr_svd.h

@@ -14,8 +14,17 @@
 
 TC_NAMESPACE_BEGIN
 
-void imp_svd(Matrix3 m, Matrix3 &u, Matrix3 &s, Matrix3 &v);
 
 void imp_svd(Matrix2 m, Matrix2 &u, Matrix2 &s, Matrix2 &v);
 
+void imp_svd(Matrix3 m, Matrix3 &u, Matrix3 &s, Matrix3 &v);
+
+void svd(Matrix2 m, Matrix2 &u, Matrix2 &sig, Matrix2 &v);
+
+void svd(Matrix3 m, Matrix3 &u, Matrix3 &sig, Matrix3 &v);
+
+void polar_decomp(Matrix2 A, Matrix2 &r, Matrix2 &s);
+
+void polar_decomp(Matrix3 A, Matrix3 &r, Matrix3 &s);
+
 TC_NAMESPACE_END

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

@@ -5,10 +5,13 @@
 if __name__ == '__main__':
     downsample = 3
     resolution = (511 / downsample, 127 / downsample, 255 / downsample)
-    tex = Texture('image', filename=tc.get_asset_path('/textures/taichi_words.png')) * 8
+
+    mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), delta_t=0.001, num_threads=8)
+
+    tex = Texture('image', filename=tc.get_asset_path('textures/taichi_words.png')) * 8
     tex = Texture('bound', tex=tex, axis=2, bounds=(0.475, 0.525), outside_val=(0, 0, 0))
-    mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), initial_velocity=(0, 0, 0), delta_t=0.001, num_threads=8,
-               density_tex=tex.id)
+    mpm.add_particles(density_tex=tex.id, initial_velocity=(0, 0, 0))
+
     for i in range(1000):
         mpm.step(0.05)
     mpm.make_video()