5 files changed, 27 insertions(+), 40 deletions(-)

M arap-mold.cpp
M arap-mold.h
M main.cpp
M mesh.cpp
M mesh.h

@@ 1,5 1,4 @@
 #include "igl/arap_rhs.h"
-#include "igl/colon.h"
 #include "igl/columnize.h"
 #include "igl/cotmatrix.h"
 #include "igl/covariance_scatter_matrix.h"


@@ 146,7 145,6 @@ bool arap::precomputation(const Mesh &mesh, arap::Data &data) {
   assert(G_sum_dim.cols() == data.CSM.rows());
   // For us, this is a no op
   data.CSM = (G_sum_dim * data.CSM).eval();
-  SDbg(data.CSM);
 
   arap_rhs(V, F, data.dim, eff_energy, data.K);
   if (flat) {


@@ 160,14 158,9 @@ bool arap::precomputation(const Mesh &mesh, arap::Data &data) {
     const double h = data.h;
     assert(h != 0);
     SparseMatrix<double> M;
-    Dbg(V);
-    Dbg(F);
     igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, data.M);
-    SDbg(data.M);
     const double dw = (1. / data.ym) * (h * h);
     SparseMatrix<double> DQ = dw * 1. / (h * h) * data.M;
-    SDbg(Q);
-    SDbg(DQ);
     Q += DQ;
     // Dummy external forces
     data.f_ext = MatrixXd::Zero(n, data.dim);


@@ 177,16 170,8 @@ bool arap::precomputation(const Mesh &mesh, arap::Data &data) {
   return min_quad_with_fixed_precompute(Q, b, SparseMatrix<double>(), true, data.solver_data);
 }
 
-// # Note: if ran with `use_dynamics` you must set `ARAPData::U0` to be your
-// initial `U` before calling this. Otherwise it'll be 0, an you will have a bad
-// time.
-void arap::iter(
-    const Mesh &mesh,
-    const arap::Args &args,
-    const Eigen::MatrixXd &bc,
-    arap::Data &data,
-    Eigen::MatrixXd &U
-    ) {
+void arap::iter(const Mesh &mesh, const arap::Args &args, const Eigen::MatrixXd &bc,
+                arap::Data &data, Eigen::MatrixXd &U) {
   using namespace Eigen;
   using namespace std;
 


@@ 227,13 212,11 @@ void arap::iter(
   // THIS NORMALIZATION IS IMPORTANT TO GET SINGLE PRECISION SVD CODE TO WORK
   // CORRECTLY.
   StackCovMats /= StackCovMats.array().abs().maxCoeff();
-  Dbg(StackCovMats);
 
   const int Rdim = data.dim;
   // This is a big matrix, which contains stacked 3x3 rotation matrices for each
   // *tet* in the mesh.
   MatrixXd Rotations(Rdim, data.CSM.rows());
-  Dbg(Rotations);
 
   if (Rdim == 2) {
     igl::fit_rotations_planar(StackCovMats, Rotations);


@@ 241,6 224,8 @@ void arap::iter(
     igl::fit_rotations(StackCovMats, true, Rotations);
   }
 
+  NO_NANS(Rotations);
+
   if (!args.skip_mold_rotations) { /* ~~ Moldability ~~ */
     // 1. Go through all the boundary faces of the tet mesh
     // 2. Skip if the face satisfies direction constraint


@@ 305,7 290,11 @@ void arap::iter(
         auto remove_dir = *data.mht_plane_n;
         if (data.orientation->at(face_i) == Orientation::Below)
           remove_dir = -remove_dir;
+        if (face_normal == -remove_dir)
+          continue;
+
         Eigen::Vector3d rotation_axis = face_normal.cross(remove_dir).normalized();
+
         auto x = remove_dir.dot(face_normal);
         auto angle_between = atan2(1.0, x);
 


@@ 421,7 410,7 @@ void arap::iter(
     // -h*data.vel = -V0+Vm1)
     // -V0-h*data.vel = -2V0+Vm1
     const double dw = (1. / data.ym) * (h * h);
-    Dl = dw * (1. / (h * h) * data.M * (-data.U0 - h * data.vel) - data.f_ext);
+    Dl = dw * (1. / (h * h) * data.M * (-mesh.rest_state - h * data.vel) - data.f_ext);
   }
 
   VectorXd Rcol;


@@ 448,6 437,6 @@ void arap::iter(
 
   if (data.with_dynamics) {
     // Keep track of velocity for next time
-    data.vel = (U - data.U0) / data.h;
+    data.vel = (U - mesh.rest_state) / data.h;
   }
 }

@@ 97,10 97,10 @@ struct Data {
 };
 
 struct Args {
-  bool should_blend;  /* whether to blend rotations with neighbors  */
-  float blend;        /* blend factor [0, 1] */
-  bool slerp_2_rots;  /* SLERP rotations to join 2. Otherwise multiply by 0.5 */
-  bool rotate_before; /* Rotate normal before checking moldability */
+  bool should_blend;        /* whether to blend rotations with neighbors  */
+  float blend;              /* blend factor [0, 1] */
+  bool slerp_2_rots;        /* SLERP rotations to join 2. Otherwise multiply by 0.5 */
+  bool rotate_before;       /* Rotate normal before checking moldability */
   bool skip_mold_rotations; /* Don't do the rotation stuff to it make moldable */
 
   Args() : should_blend(false), blend(0.0), slerp_2_rots(false), rotate_before(false){};


@@ 116,19 116,13 @@ struct Args {
 //   b  #b list of "boundary" fixed vertex indices into V
 // Outputs:
 //   data  struct containing necessary precomputation
-bool precomputation(const Mesh &, Data &data); 
+bool precomputation(const Mesh &, Data &data);
 
 // Inputs:
 //   bc  #b by dim list of boundary conditions
 //   data  struct containing necessary precomputation and parameters
 //   U  #V by dim initial guess
-void iter(
-    const Mesh &,
-    const Args &,
-    const Eigen::MatrixXd &bc,
-    Data &data,
-    Eigen::MatrixXd &U
-    );
+void iter(const Mesh &, const Args &, const Eigen::MatrixXd &bc, Data &data, Eigen::MatrixXd &U);
 
 Moldable::Enum is_moldable(Eigen::Vector3d &remove_dir, std::vector<bool> &is_point_above,
                            Eigen::MatrixXi &faces, int face, Eigen::Vector3d normal,

@@ 57,6 57,8 @@ bool show_forces = false;
 // Don't use the ARAP enforcing moldability.  (only makes sense to use wiht `stitching`).
 bool regular_arap = false;
 
+int iteration_count = 0;
+
 // Cut plane and point with `float`, for ImGui.
 Eigen::Vector3f fplane_p;
 Eigen::Vector3f fplane_n;


@@ 102,6 104,7 @@ void draw_menu(Mesh &mesh) {
     if (ImGui::Button("Stop ARAP"))
       run_arap = false;
   }
+  ImGui::Text("iter count: %d", iteration_count);
 
   if (ImGui::Button("Reset")) {
     mesh.vertices = mesh.rest_state;


@@ 172,6 175,7 @@ void draw_menu(Mesh &mesh) {
 
 // Perform one ARAP step.
 void step_arap(Mesh &mesh) {
+  iteration_count++;
   Eigen::MatrixXd bc(BOUNDARY_SIZE, 3);
   for (int i = 0; i < BOUNDARY_SIZE; i++)
     bc.row(i) = mesh.rest_state.row(mesh.boundary(i));


@@ 210,7 214,8 @@ void compute_things(Mesh &mesh, bool only_plane) {
     for (int f = 0; f < mesh.faces.rows(); f++) {
       double dot;
       Eigen::Vector3d removal_dir = plane_n;
-      moldable[f] = arap::is_moldable(removal_dir, is_point_above, mesh.faces, f, mesh.normals.row(f), &dot);
+      moldable[f] =
+          arap::is_moldable(removal_dir, is_point_above, mesh.faces, f, mesh.normals.row(f), &dot);
       if (moldable[f] == arap::Moldable::Crossing) {
         orientation[f] = arap::Orientation::Crossing;
       } else {


@@ 353,7 358,7 @@ void draw_plane(Viewer &viewer, Mesh &mesh) {
 
   if (show_forces) {
     viewer.data().add_edges(mesh.vertices, mesh.vertices + mesh.forces,
-        Eigen::RowVector3d(1.0, 0.0, 1.0));
+                            Eigen::RowVector3d(1.0, 0.0, 1.0));
   }
 }
 


@@ 369,7 374,8 @@ int ray_intersect(Ray ray, int src_vertex, Mesh &mesh) {
   int best_vertex = -1;
 
   for (int f = 0; f < mesh.faces.rows(); f++) {
-    int res = intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.faces, &collision_point);
+    int res = intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.faces,
+                                      &collision_point);
     if (res != 1)
       continue;
 


@@ 505,8 511,8 @@ int main(int argc, char *argv[]) {
   arap_data.mht_tf2t = mesh.face2tets;
   arap_data.mht_plane_p = &plane_p;
   arap_data.mht_plane_n = &plane_n;
-  arap_data.mht_tet_rotations =
-      std::vector<Eigen::Quaternion<double>>(mesh.tets.rows(), Eigen::Quaternion<double>(1, 0, 0, 0));
+  arap_data.mht_tet_rotations = std::vector<Eigen::Quaternion<double>>(
+      mesh.tets.rows(), Eigen::Quaternion<double>(1, 0, 0, 0));
   arap_data.vert_is_above = &is_point_above;
   arap_data.moldable = &moldable;
   arap_data.orientation = &orientation;

@@ 172,4 172,3 @@ void mesh_print_sizes(Mesh &mesh) {
   Dbg(mesh.boundary);
   Dbg(mesh.forces);
 }
-

@@ 3,7 3,6 @@
 #define Dbg(x) printf("%20s: %6zu by %6zu\n", #x, x.rows(), x.cols())
 #define SDbg(x) printf("%20s: %6d by %6d\n", #x, x.rows(), x.cols())
 
-
 #include <Eigen/Core>
 
 using Eigen::MatrixXd;