5 files changed, 201 insertions(+), 195 deletions(-)

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

@@ 54,11 54,15 @@ arap::Moldable::Enum arap::is_moldable(Eigen::Vector3d &remove_dir,
   return Moldable::Yes;
 }
 
-bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
-                          const Eigen::MatrixXi &b, arap::Data &data) {
+bool arap::precomputation(const Mesh &mesh, arap::Data &data) {
   using namespace std;
   using namespace Eigen;
-  typedef typename DerivedV::Scalar Scalar;
+
+  const Eigen::MatrixXd &V = mesh.vertices;
+  const Eigen::MatrixXi &F = mesh.tets;
+  const int dim = mesh.vertices.cols();
+  const Eigen::VectorXi &b = mesh.boundary;
+
   // number of vertices
   const int n = V.rows();
   data.n = n;


@@ 75,19 79,16 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
   assert(data.dim <= V.cols() && "solve dim should be <= embedding");
   bool flat = (V.cols() - data.dim) == 1;
 
-  DerivedV plane_V;
-  DerivedF plane_F;
-  typedef SparseMatrix<Scalar> SparseMatrixS;
+  MatrixXd plane_V;
+  MatrixXi plane_F;
 
-  SparseMatrixS ref_map, ref_map_dim;
+  SparseMatrix<double> ref_map, ref_map_dim;
   if (flat) {
     igl::project_isometrically_to_plane(V, F, plane_V, plane_F, ref_map);
     igl::repdiag(ref_map, dim, ref_map_dim);
   }
 
-  const PlainObjectBase<DerivedV> &ref_V = (flat ? plane_V : V);
-  const PlainObjectBase<DerivedF> &ref_F = (flat ? plane_F : F);
-  SparseMatrixS L;
+  SparseMatrix<double> L;
   igl::cotmatrix(V, F, L);
 
   igl::ARAPEnergyType eff_energy = data.energy;


@@ 110,7 111,7 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
 
   // Get covariance scatter matrix, when applied collects the covariance
   // matrices used to fit rotations to during optimization
-  covariance_scatter_matrix(ref_V, ref_F, eff_energy, data.CSM);
+  covariance_scatter_matrix(V, F, eff_energy, data.CSM);
   if (flat) {
     data.CSM = (data.CSM * ref_map_dim.transpose()).eval();
   }


@@ 145,8 146,9 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
   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(ref_V, ref_F, data.dim, eff_energy, data.K);
+  arap_rhs(V, F, data.dim, eff_energy, data.K);
   if (flat) {
     data.K = (ref_map_dim * data.K).eval();
   }


@@ 158,9 160,14 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
     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);


@@ 173,10 180,19 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
 // # 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 arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &data, int num_tets,
-                Eigen::MatrixXi &TNeigh, 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;
+
+  const MatrixXi &TNeigh = mesh.neighbors;
+  const int num_tets = mesh.tets.rows();
+
   assert(data.b.size() == bc.rows());
   if (bc.size() > 0) {
     assert(bc.cols() == data.dim && "bc.cols() match data.dim");


@@ 211,11 227,13 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
   // 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);

@@ 5,6 5,8 @@
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 
+#include "mesh.h"
+
 namespace arap {
 
 typedef Eigen::MatrixXd DerivedV;


@@ 114,15 116,19 @@ struct Args {
 //   b  #b list of "boundary" fixed vertex indices into V
 // Outputs:
 //   data  struct containing necessary precomputation
-bool precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
-                    const Eigen::MatrixXi &b, 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 Args &, const Eigen::MatrixXd &bc, Data &data, int num_tets, Eigen::MatrixXi &,
-          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,

@@ 12,7 12,6 @@
 #include "mesh.h"
 
 #define BOUNDARY_SIZE 4
-#define Dbg(x) printf("%10s: %7zu by %3zu\n", #x, x.rows(), x.cols())
 
 using igl::opengl::glfw::Viewer;
 


@@ 26,41 25,14 @@ struct Ray {
   Vector3d dir;
 };
 
-void step_arap();
-bool pre_draw(Viewer &);
-void draw_plane(Viewer &);
-bool key_down(Viewer &, unsigned char, int);
-void draw_menu();
+void step_arap(Mesh &);
+void draw_plane(Viewer &, Mesh &);
+void draw_menu(Mesh &);
 void update_vectors();
-void ensure_moldability();
-void update_rest_state();
-void remesh();
-int ray_intersect(Ray, int, MatrixXd &, MatrixXi &);
-void compute_things(MatrixXd &, MatrixXi &, bool = false);
+int ray_intersect(Ray, int, Mesh &);
+void compute_things(Mesh &, bool = false);
 int intersect3D_RayTriangle(Ray, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
 
-// Model data
-MatrixXd V;
-MatrixXi F;
-MatrixXd colors;
-
-// Tetrahedralized data
-MatrixXd TV;     // Vertices
-MatrixXi TT;     // Tets
-MatrixXi TF;     // Triangles
-MatrixXd TN;     // Normals
-MatrixXi TF2T;   // Face to tets mapping
-MatrixXi TT2F;   // Tet to faces mapping
-MatrixXi TNeigh; // Neighbor tets for each tet
-
-MatrixXd RestState;
-// Indices of the fixed vertices in TV.
-VectorXi boundary(BOUNDARY_SIZE);
-// ARAP'd vertices
-MatrixXd U;  /* Deformed vertices */
-MatrixXd U0; /* Original U  (this is really just TV) */
-MatrixXd Forces;
-
 // Various computed things about the mesh
 // If a point is above the plane.
 std::vector<bool> is_point_above;


@@ 105,13 77,8 @@ void update_vectors() {
   plane_n = fplane_n.cast<double>().normalized();
 }
 
-void update_rest_state() {
-  RestState = U;
-  arap::precomputation(RestState, TT, RestState.cols(), boundary, arap_data);
-}
-
 // Draw the custom menu
-void draw_menu() {
+void draw_menu(Mesh &mesh) {
   bool recompute = false;
   bool px = ImGui::SliderFloat("plane_p.x", &fplane_p[0], mins.x(), maxs.x(), "%.4f");
   bool py = ImGui::SliderFloat("plane_p.y", &fplane_p[1], mins.y(), maxs.y(), "%.4f");


@@ 137,17 104,18 @@ void draw_menu() {
   }
 
   if (ImGui::Button("Reset")) {
-    U = RestState;
+    mesh.vertices = mesh.rest_state;
     recompute |= true;
   }
 
   if (ImGui::Button("Update Rest State")) {
-    update_rest_state();
+    mesh.rest_state = mesh.vertices;
+    arap::precomputation(mesh, arap_data);
     recompute |= true;
   }
 
   if (ImGui::Button("Step ARAP")) {
-    step_arap();
+    step_arap(mesh);
   }
 
   ImGui::Checkbox("Show Forces", &show_forces);


@@ 159,7 127,7 @@ void draw_menu() {
 
   if (ImGui::Checkbox("Stiching", &stitching)) {
     if (!stitching)
-      Forces = MatrixXd::Zero(U.rows(), 3);
+      mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), mesh.vertices.cols());
     recompute |= true;
   }
   ImGui::Checkbox("Skip Moldability Rotations", &regular_arap);


@@ 196,33 164,21 @@ void draw_menu() {
     }
     ImGui::EndChild();
   }
-  if (recompute) {
-    compute_things(U, TF);
-  }
-}
 
-// Handle key down events.
-bool key_down(Viewer &viewer, unsigned char key, int mod) {
-  if (key == ' ') {
-    run_arap = !run_arap;
-  } else if (key == '.') {
-    step_arap();
-  } else if (key == 'R') {
-    RestState = U = TV;
-    compute_things(U, TF);
+  if (recompute) {
+    compute_things(mesh);
   }
-  return false;
 }
 
 // Perform one ARAP step.
-void step_arap() {
+void step_arap(Mesh &mesh) {
   Eigen::MatrixXd bc(BOUNDARY_SIZE, 3);
   for (int i = 0; i < BOUNDARY_SIZE; i++)
-    bc.row(i) = RestState.row(boundary(i));
+    bc.row(i) = mesh.rest_state.row(mesh.boundary(i));
 
   arap_args.skip_mold_rotations = regular_arap;
-  arap::iter(arap_args, bc, arap_data, TT.rows(), TNeigh, U);
-  compute_things(U, TF);
+  arap::iter(mesh, arap_args, bc, arap_data, mesh.vertices);
+  compute_things(mesh);
 }
 
 // Compute all the things that we need each time the vertex set changes. Here's


@@ 230,31 186,31 @@ void step_arap() {
 //  - `TN` (normals)
 //  - `is_point_above`
 //  - `colors`
-void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only_plane) {
+void compute_things(Mesh &mesh, bool only_plane) {
 
   if (!only_plane) { /* `normals`: depends on vertices */
-    igl::per_face_normals(vertices, faces, TN);
-    viewer.data().set_vertices(vertices);
-    viewer.data().set_normals(TN);
+    igl::per_face_normals(mesh.vertices, mesh.faces, mesh.normals);
+    viewer.data().set_vertices(mesh.vertices);
+    viewer.data().set_normals(mesh.normals);
   }
 
   { /* `is_point_above`: depends on vertices and the plane */
-    is_point_above.resize(vertices.rows());
-    for (int v_i = 0; v_i < vertices.rows(); v_i++) {
-      Vector3d v(vertices.row(v_i));
+    is_point_above.resize(mesh.vertices.rows());
+    for (int v_i = 0; v_i < mesh.vertices.rows(); v_i++) {
+      Vector3d v(mesh.vertices.row(v_i));
       is_point_above[v_i] = (v - plane_p).dot(plane_n) >= 0.0;
     }
   }
 
   { /* `moldable`, `orientation`, and `dots`: depends on vertices and the plane
      */
-    moldable.resize(faces.rows());
-    orientation.resize(faces.rows());
-    dots.resize(faces.rows());
-    for (int f = 0; f < faces.rows(); f++) {
+    moldable.resize(mesh.faces.rows());
+    orientation.resize(mesh.faces.rows());
+    dots.resize(mesh.faces.rows());
+    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, faces, f, TN.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 {


@@ 269,38 225,38 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
   }
 
   { /* `colors`: depends on dot */
-    colors.resize(faces.rows(), 3);
-    for (int f = 0; f < faces.rows(); f++) {
+    mesh.face_colors.resize(mesh.faces.rows(), 3);
+    for (int f = 0; f < mesh.faces.rows(); f++) {
       switch (moldable[f]) {
       case arap::Moldable::Yes: {
-        colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
+        mesh.face_colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
       } break;
       case arap::Moldable::Crossing: {
-        colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
+        mesh.face_colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
       } break;
       default: {
         double r = -dots[f];
         r = sqrt(r);
-        colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
+        mesh.face_colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
       } break;
       }
     }
-    viewer.data().set_colors(colors);
+    viewer.data().set_colors(mesh.face_colors);
   }
 
   if (stitching) {
-    Forces = MatrixXd::Zero(vertices.rows(), 3);
-    std::vector<bool> seen(vertices.rows());
-    for (int f = 0; f < faces.rows(); f++) {
+    mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), 3);
+    std::vector<bool> seen(mesh.vertices.rows());
+    for (int f = 0; f < mesh.faces.rows(); f++) {
       auto m = moldable[f];
       // XXX(todo): what do we do about `Crossing`? For now, skip as well.
       if (m == arap::Moldable::Yes)
         continue;
       if (m == arap::Moldable::Crossing) {
         for (int i = 0; i < 3; i++) {
-          int v = faces(f, i);
+          int v = mesh.faces(f, i);
           seen[v] = true;
-          Forces.row(v) = Vector3d::Zero();
+          mesh.forces.row(v) = Vector3d::Zero();
         }
         continue;
       }


@@ 308,50 264,39 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
           (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
 
       for (int i = 0; i < 3; i++) {
-        int v = faces(f, i);
+        int v = mesh.faces(f, i);
         // Only do a vertex one time; this is okay since the direction is +- plane_n, and not the
         // triangle normal.
         if (seen[v])
           continue;
         seen[v] = true;
 
-        Vector3d vertex = vertices.row(v);
+        Vector3d vertex = mesh.vertices.row(v);
         Ray ray;
         ray.src = vertex;
         ray.dir = ray_dir;
 
-        int other_i = ray_intersect(ray, v, vertices, faces);
+        int other_i = ray_intersect(ray, v, mesh);
         if (other_i == -1)
           continue;
 
-        Vector3d other_vertex = vertices.row(other_i);
+        Vector3d other_vertex = mesh.vertices.row(other_i);
         auto force_direction = other_vertex - vertex;
         // XXX(todo): what should the magnitude of the force be? Should probably take a closer look
         // at the dynamics implementation in `arap-mold.cpp`.
 
         // Symmetric forces - think springs
-        Forces.row(v) += force_direction * 0.48;
-        Forces.row(other_i) -= force_direction * 0.48;
+        mesh.forces.row(v) += force_direction * 0.48;
+        mesh.forces.row(other_i) -= force_direction * 0.48;
       }
     }
 
-    arap_data.f_ext = Forces;
+    arap_data.f_ext = mesh.forces;
   }
 }
 
-// This is called before each draw call. Compute normals and set right colors.
-// XXX(note): much of this should probably just be done after `step_arap`, since
-// the rest is pretty constant.
-bool pre_draw(Viewer &viewer) {
-  draw_plane(viewer);
-  if (run_arap) {
-    step_arap();
-  }
-  return false;
-}
-
 // Draw the cut plane in a very hacky way.
-void draw_plane(Viewer &viewer) {
+void draw_plane(Viewer &viewer, Mesh &mesh) {
   // Must clear the previous plane from the screen. Maybe there's a better way
   // of doing this than to clear the entire thing.
   viewer.data().lines.resize(0, 0);


@@ 407,14 352,15 @@ void draw_plane(Viewer &viewer) {
   viewer.data().add_edges(planeEA, planeEB, c);
 
   if (show_forces) {
-    viewer.data().add_edges(U, U + Forces, Eigen::RowVector3d(1.0, 0.0, 1.0));
+    viewer.data().add_edges(mesh.vertices, mesh.vertices + mesh.forces,
+        Eigen::RowVector3d(1.0, 0.0, 1.0));
   }
 }
 
 // Shoot a ray from `p` in direction `v`, and see if it intersects with any of
 // the `faces`. If it does, return the index of the closest vertex to the
 // intersection point. If not, return `-1`.
-int ray_intersect(Ray ray, int src_vertex, Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces) {
+int ray_intersect(Ray ray, int src_vertex, Mesh &mesh) {
   // XXX(todo): Option to intersect with cut plane. Maybe this isn't a good idea - rays from any
   //            non-moldable face will intersect with the plane!
 


@@ 422,8 368,8 @@ int ray_intersect(Ray ray, int src_vertex, Eigen::MatrixXd &vertices, Eigen::Mat
   double best_dist = std::numeric_limits<double>::max();
   int best_vertex = -1;
 
-  for (int f = 0; f < faces.rows(); f++) {
-    int res = intersect3D_RayTriangle(ray, TN.row(f), f, vertices, faces, &collision_point);
+  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);
     if (res != 1)
       continue;
 


@@ 432,21 378,21 @@ int ray_intersect(Ray ray, int src_vertex, Eigen::MatrixXd &vertices, Eigen::Mat
 
     bool is_own_face = false;
     for (int i = 0; i < 3; i++)
-      if (faces(f, i) == src_vertex)
+      if (mesh.faces(f, i) == src_vertex)
         is_own_face = true;
     if (is_own_face)
       continue;
 
     for (int i = 0; i < 3; i++) {
-      int v = faces(f, i);
-      double norm = (collision_point - Vector3d(vertices.row(v))).norm();
+      int v = mesh.faces(f, i);
+      double norm = (collision_point - Vector3d(mesh.vertices.row(v))).norm();
       if (norm < dist_to_collision) {
         dist_to_collision = norm;
         closest_vertex = v;
       }
     }
 
-    double dist_to_ray_pos = (Vector3d(vertices.row(closest_vertex)) - ray.src).norm();
+    double dist_to_ray_pos = (Vector3d(mesh.vertices.row(closest_vertex)) - ray.src).norm();
     if (dist_to_ray_pos < best_dist) {
       best_dist = dist_to_ray_pos;
       best_vertex = closest_vertex;


@@ 512,28 458,6 @@ int intersect3D_RayTriangle(Ray ray, Vector3d triangle_normal, int triangle, Mat
   return 1; // I is in T
 }
 
-// // Remesh the model based on the deformed version.
-// // WARNING: THIS WILL DESTROY THE ORIGINAL MODEL, SO `reset` DOES NOTHING!
-// void remesh() {  NOTE: This is not remeshing, but triangulating..
-//   TetmeshIn tetmesh_in;
-//   TetmeshOut tetmesh_out;
-//   tetmesh_in.vertices = &U;
-//   tetmesh_in.faces = &TF;
-//
-//   int error = make_tetmesh(tetmesh_in, tetmesh_out);
-//   if (error) {
-//     std::cerr << "[ERR] tetesh returned " << error << std::endl;
-//     exit(1);
-//   }
-//
-//   TV = tetmesh_out.vertices;
-//   U = TV;
-//   TT = tetmesh_out.tets;
-//   TF = tetmesh_out.faces;
-//   TF2T = tetmesh_out.face2tets;
-//   TT2F = tetmesh_out.tet2faces;
-// }
-
 int main(int argc, char *argv[]) {
   int error;
   using namespace std;


@@ 543,57 467,57 @@ int main(int argc, char *argv[]) {
     return 1;
   }
 
-  igl::readOBJ(argv[1], V, F);
-
-  TetmeshIn tetmesh_in;
   TetmeshOut tetmesh_out;
-  tetmesh_in.vertices = &V;
-  tetmesh_in.faces = &F;
-
-  error = make_tetmesh(tetmesh_in, tetmesh_out);
-  if (error) {
-    cerr << "[ERR] tetesh returned " << error << endl;
-    exit(1);
+  {
+    MatrixXd read_verts;
+    MatrixXi read_faces;
+    igl::readOBJ(argv[1], read_verts, read_faces);
+
+    TetmeshIn tetmesh_in;
+    tetmesh_in.vertices = &read_verts;
+    tetmesh_in.faces = &read_faces;
+
+    error = make_tetmesh(tetmesh_in, tetmesh_out);
+    if (error) {
+      cerr << "[ERR] tetesh returned " << error << endl;
+      exit(1);
+    }
   }
 
-  TV = tetmesh_out.vertices;
-  TT = tetmesh_out.tets;
-  TF = tetmesh_out.faces;
-  TF2T = tetmesh_out.face2tets;
-  TT2F = tetmesh_out.tet2faces;
-  TNeigh = tetmesh_out.neighbors;
-  RestState = U = TV;
-  Forces = MatrixXd::Zero(TV.rows(), 3);
-
-  boundary << 0, 1, 2, 3; // 19, 22, 121, 122, 221, 224, 321, 326;
-
-  //   Dbg(V);
-  //   Dbg(F);
-  //   Dbg(TV);
-  //   Dbg(TT);
-  //   Dbg(TF);
-  //   Dbg(TT2F);
-  //   Dbg(TF2T);
-  //   Dbg(TNeigh);
-
-  arap_data.faces = &TF;
-  arap_data.normals = &TN;
-  arap_data.mht_tf2t = TF2T;
+  Mesh mesh;
+
+  mesh.vertices = tetmesh_out.vertices;
+  mesh.faces = tetmesh_out.faces;
+  mesh.tets = tetmesh_out.tets;
+  mesh.face2tets = tetmesh_out.face2tets;
+  mesh.tet2faces = tetmesh_out.tet2faces;
+  mesh.neighbors = tetmesh_out.neighbors;
+
+  mesh.rest_state = mesh.vertices;
+  mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), mesh.vertices.cols());
+  mesh.boundary = VectorXi(4);
+  mesh.boundary << 0, 1, 2, 3;
+
+  mesh_print_sizes(mesh);
+
+  arap_data.faces = &mesh.faces;
+  arap_data.normals = &mesh.normals;
+  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>>(TT.size(), Eigen::Quaternion<double>(1, 0, 0, 0));
+      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;
-  arap_data.U0 = TV;
+  arap_data.U0 = mesh.vertices;
   arap_data.with_dynamics = true;
   // XXX(todo): what's a reasonable value here?
   arap_data.h = 0.01;
-  arap::precomputation(TV, TT, TV.cols(), boundary, arap_data);
+  arap::precomputation(mesh, arap_data);
 
-  maxs = V.colwise().maxCoeff().cast<float>();
-  mins = V.colwise().minCoeff().cast<float>();
+  maxs = mesh.vertices.colwise().maxCoeff().cast<float>();
+  mins = mesh.vertices.colwise().minCoeff().cast<float>();
   scale = (maxs - mins).minCoeff();
 
   fplane_p = ((maxs + mins) / 2.0);


@@ 603,17 527,36 @@ int main(int argc, char *argv[]) {
   igl::opengl::glfw::imgui::ImGuiMenu menu;
   viewer.plugins.push_back(&menu);
 
-  menu.callback_draw_custom_window = &draw_menu;
-  menu.callback_draw_viewer_menu = [&]() { menu.draw_viewer_menu(); };
+  menu.callback_draw_custom_window = [&mesh]() { draw_menu(mesh); };
+  menu.callback_draw_viewer_menu = [&menu]() { menu.draw_viewer_menu(); };
+
+  viewer.callback_pre_draw = [&mesh](Viewer &viewer) {
+    draw_plane(viewer, mesh);
+    if (run_arap) {
+      step_arap(mesh);
+    }
+    return false;
+  };
+
+  viewer.callback_key_down = [&mesh, &tetmesh_out](Viewer &viewer, unsigned char key, int mod) {
+    if (key == ' ') {
+      run_arap = !run_arap;
+    } else if (key == '.') {
+      step_arap(mesh);
+    } else if (key == 'R') {
+      mesh.vertices = tetmesh_out.vertices;
+      mesh.rest_state = mesh.vertices;
+      compute_things(mesh);
+    }
+    return false;
+  };
 
-  viewer.callback_pre_draw = &pre_draw;
-  viewer.callback_key_down = &key_down;
   viewer.data().clear();
-  viewer.data().set_mesh(TV, TF);
+  viewer.data().set_mesh(mesh.vertices, mesh.faces);
   viewer.data().set_face_based(true);
   viewer.core().is_animating = true;
 
-  compute_things(TV, TF);
+  compute_things(mesh);
 
   viewer.launch();
 }

@@ 158,3 158,18 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
 
   return 0;
 }
+
+void mesh_print_sizes(Mesh &mesh) {
+  Dbg(mesh.vertices);
+  Dbg(mesh.faces);
+  Dbg(mesh.tets);
+  Dbg(mesh.normals);
+  Dbg(mesh.face2tets);
+  Dbg(mesh.tet2faces);
+  Dbg(mesh.neighbors);
+  Dbg(mesh.face_colors);
+  Dbg(mesh.rest_state);
+  Dbg(mesh.boundary);
+  Dbg(mesh.forces);
+}
+

@@ 1,9 1,14 @@
 #pragma once
 
+#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;
 using Eigen::MatrixXi;
+using Eigen::VectorXi;
 
 // Input here is a triangle mesh.
 struct TetmeshIn {


@@ 15,6 20,7 @@ struct TetmeshIn {
   TetmeshIn(){};
 };
 
+// Output from tetmesh
 struct TetmeshOut {
   // Vertices
   MatrixXd vertices;


@@ 32,6 38,24 @@ struct TetmeshOut {
   TetmeshOut(){};
 };
 
+// All data for our mesh
+struct Mesh {
+  MatrixXd vertices;
+  MatrixXi faces;
+  MatrixXi tets;
+  MatrixXd normals;
+  MatrixXi face2tets;
+  MatrixXi tet2faces;
+  MatrixXi neighbors;
+
+  MatrixXd face_colors;
+  MatrixXd rest_state;
+  VectorXi boundary;
+  MatrixXd forces;
+};
+
+void mesh_print_sizes(Mesh &);
+
 // Turns a 3D triangle mesh into a tetmesh, using `tetgen`.
 // Returns 0 if OK.
 int make_tetmesh(TetmeshIn &in, TetmeshOut &out);