dynamics in arap

Former-commit-id: c0bb7ee84e85b7afc7a36edf533a80f0f7d2fa76

examples/arap/example.cpp

@@ -69,6 +69,9 @@ Eigen::Quaterniond animation_to_quat;
 // Use vector for range-based `for`
 std::vector<State> undo_stack;
 std::vector<State> redo_stack;
+bool paused = false;
+double t = 0;
+double pause_time = 0.0;
 
 void push_undo()
 {
@@ -211,7 +214,10 @@ bool init_arap()
   {
     return false;
   }
-  partition(W,100,arap_data.G,_S,_D);
+  //arap_data.with_dynamics = true;
+  //arap_data.h = 0.5;
+  //arap_data.max_iter = 100;
+  //partition(W,100,arap_data.G,_S,_D);
   return arap_precomputation(V,*Ele,b,arap_data);
 }
 
@@ -223,6 +229,10 @@ bool update_arap()
   MatrixXd bc(num_in_selection(S),V.cols());
   // get b from S
   {
+    if(!paused)
+    {
+      t = get_seconds()-pause_time;
+    }
     int bi = 0;
     for(int v = 0;v<S.rows(); v++)
     {
@@ -234,15 +244,15 @@ bool update_arap()
           case 0:
           {
             const double r = mid(0)*0.25;
-            bc(bi,0) += r*cos(0.5*get_seconds()*2.*PI);
-            bc(bi,1) -= r+r*sin(0.5*get_seconds()*2.*PI);
+            bc(bi,0) += r*cos(0.5*t*2.*PI);
+            bc(bi,1) -= r+r*sin(0.5*t*2.*PI);
             break;
           }
           case 1:
           {
             const double r = mid(1)*0.15;
-            bc(bi,1) += r+r*cos(0.15*get_seconds()*2.*PI);
-            bc(bi,2) -= r*sin(0.15*get_seconds()*2.*PI);
+            bc(bi,1) += r+r*cos(0.15*t*2.*PI);
+            bc(bi,2) -= r*sin(0.15*t*2.*PI);
 
             //// Pull-up
             //bc(bi,0) += 0.42;//mid(0)*0.5;
@@ -586,6 +596,7 @@ void mouse_drag(int mouse_x, int mouse_y)
 void key(unsigned char key, int mouse_x, int mouse_y)
 {
   using namespace std;
+  using namespace igl;
   switch(key)
   {
     // ESC
@@ -594,6 +605,20 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     // ^C
     case char(3):
       exit(0);
+    case ' ':
+      {
+        static double pause_start,pause_stop;
+        paused = !paused;
+        if(paused)
+        {
+          pause_start = get_seconds();
+        }else
+        {
+          pause_stop = get_seconds();
+          pause_time += (pause_stop-pause_start);
+        }
+        break;
+      }
     default:
       if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
       {

include/igl/massmatrix.cpp

@@ -21,12 +21,21 @@ IGL_INLINE void igl::massmatrix(
 {
   using namespace Eigen;
   using namespace std;
-  assert(type!=MASSMATRIX_FULL);
 
   const int n = V.rows();
   const int m = F.rows();
   const int simplex_size = F.cols();
 
+  MassMatrixType eff_type = type;
+  // Use voronoi of for triangles by default, otherwise barycentric
+  if(type == MASSMATRIX_DEFAULT)
+  {
+    eff_type = (simplex_size == 3?MASSMATRIX_VORONOI:MASSMATRIX_BARYCENTRIC);
+  }
+
+  // Not yet supported
+  assert(type!=MASSMATRIX_FULL);
+
   Matrix<int,Dynamic,1> MI;
   Matrix<int,Dynamic,1> MJ;
   Matrix<Scalar,Dynamic,1> MV;
@@ -52,7 +61,7 @@ IGL_INLINE void igl::massmatrix(
       dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
     }
 
-    switch(type)
+    switch(eff_type)
     {
       case MASSMATRIX_BARYCENTRIC:
         // diagonal entries for each face corner
@@ -115,13 +124,13 @@ IGL_INLINE void igl::massmatrix(
         assert(false && "Implementation incomplete");
         break;
       default:
-        assert(false && "Unknown Mass matrix type");
+        assert(false && "Unknown Mass matrix eff_type");
     }
 
   }else if(simplex_size == 4)
   {
     assert(V.cols() == 3);
-    assert(type == MASSMATRIX_BARYCENTRIC);
+    assert(eff_type == MASSMATRIX_BARYCENTRIC);
     MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
     MI.block(0*m,0,m,1) = F.col(0);
     MI.block(1*m,0,m,1) = F.col(1);

include/igl/massmatrix.h

@@ -18,11 +18,12 @@ namespace igl
 
   enum MassMatrixType
   {
-    MASSMATRIX_BARYCENTRIC,
-    MASSMATRIX_VORONOI,
-    MASSMATRIX_FULL
+    MASSMATRIX_BARYCENTRIC = 0,
+    MASSMATRIX_VORONOI = 1,
+    MASSMATRIX_FULL = 2,
+    MASSMATRIX_DEFAULT = 3,
+    NUM_MASSMATRIX = 4
   };
-#define NUM_MASSMATRIXTYPE 3
 
   // Constructs the mass (area) matrix for a given mesh (V,F).
   //

include/igl/svd3x3/arap.cpp

@@ -8,6 +8,7 @@
 #include "arap.h"
 #include <igl/colon.h>
 #include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
 #include <igl/group_sum_matrix.h>
 #include <igl/covariance_scatter_matrix.h>
 #include <igl/speye.h>
@@ -113,9 +114,24 @@ IGL_INLINE bool igl::arap_precomputation(
   assert(G_sum_dim.cols() == data.CSM.rows());
   data.CSM = (G_sum_dim * data.CSM).eval();
 
+
   arap_rhs(V,F,eff_energy,data.K);
 
   SparseMatrix<double> Q = (-0.5*L).eval();
+
+  if(data.with_dynamics)
+  {
+    const double h = data.h;
+    assert(h != 0);
+    SparseMatrix<double> M;
+    massmatrix(V,F,MASSMATRIX_DEFAULT,data.M);
+    SparseMatrix<double> DQ = 0.5/(h*h)*data.M;
+    Q += DQ;
+    // Dummy external forces
+    data.f_ext = MatrixXd::Zero(n,dim);
+    data.vel = MatrixXd::Zero(n,dim);
+  }
+
   return min_quad_with_fixed_precompute(
     Q,b,SparseMatrix<double>(),true,data.solver_data);
 }
@@ -140,7 +156,14 @@ IGL_INLINE bool igl::arap_solve(
     // terrible initial guess.. should at least copy input mesh
     U = MatrixXd::Zero(data.n,dim);
   }
+  // changes each arap iteration
   MatrixXd U_prev = U;
+  // doesn't change for fixed with_dynamics timestep
+  MatrixXd U0;
+  if(data.with_dynamics)
+  {
+    U0 = U_prev;
+  }
   while(iter < data.max_iter)
   {
     U_prev = U;
@@ -181,6 +204,21 @@ IGL_INLINE bool igl::arap_solve(
       }
     }
 
+    MatrixXd Dl;
+    if(data.with_dynamics)
+    {
+      assert(M.rows() == n && 
+        "No mass matrix. Call arap_precomputation if changing with_dynamics");
+      const double h = data.h;
+      assert(h != 0);
+      //Dl = 1./(h*h*h)*M*(-2.*V0 + Vm1) - fext;
+      // data.vel = (V0-Vm1)/h
+      // h*data.vel = (V0-Vm1)
+      // -h*data.vel = -V0+Vm1)
+      // -V0-h*data.vel = -2V0+Vm1
+      Dl = 1./(h*h)*data.M*(-U0 - h*data.vel) - data.f_ext;
+    }
+
     VectorXd Rcol;
     columnize(eff_R,num_rots,2,Rcol);
     VectorXd Bcol = -data.K * Rcol;
@@ -188,6 +226,10 @@ IGL_INLINE bool igl::arap_solve(
     {
       VectorXd Uc,Bc,bcc,Beq;
       Bc = Bcol.block(c*n,0,n,1);
+      if(data.with_dynamics)
+      {
+        Bc += Dl.col(c);
+      }
       bcc = bc.col(c);
       min_quad_with_fixed_solve(
         data.solver_data,
@@ -195,10 +237,14 @@ IGL_INLINE bool igl::arap_solve(
         Uc);
       U.col(c) = Uc;
     }
-    
 
     iter++;
   }
+  if(data.with_dynamics)
+  {
+    // Keep track of velocity for next time
+    data.vel = (U-U0)/data.h;
+  }
   return true;
 }
 

include/igl/svd3x3/arap.h

@@ -21,21 +21,24 @@ namespace igl
     // G  #V list of group indices (1 to k) for each vertex, such that vertex i
     //    is assigned to group G(i)
     // energy  type of energy to use
-    // with_dynamics  whether using dynamics 
+    // with_dynamics  whether using dynamics (need to call arap_precomputation
+    //   after changing)
     // f_ext  #V by dim list of external forces
+    // vel  #V by dim list of velocities
     // h  dynamics time step
     // max_iter  maximum inner iterations
     // K  rhs pre-multiplier
+    // M  mass matrix
     // solver_data  quadratic solver data
     // b  list of boundary indices into V
     int n;
     Eigen::VectorXi G;
     ARAPEnergyType energy;
     bool with_dynamics;
-    Eigen::MatrixXd f_ext;
+    Eigen::MatrixXd f_ext,vel;
     double h;
     int max_iter;
-    Eigen::SparseMatrix<double> K;
+    Eigen::SparseMatrix<double> K,M;
     Eigen::SparseMatrix<double> CSM;
     min_quad_with_fixed_data<double> solver_data;
     Eigen::VectorXi b;