Feilong Song HW4

src/arap_precompute.cpp

@@ -1,4 +1,8 @@
 #include "arap_precompute.h"
+#include <igl/cotmatrix.h>
+#include <igl/cotmatrix_entries.h>
+#include <igl/min_quad_with_fixed.h>
+#include <vector>
 
 void arap_precompute(
   const Eigen::MatrixXd & V,
@@ -7,5 +11,33 @@ void arap_precompute(
   igl::min_quad_with_fixed_data<double> & data,
   Eigen::SparseMatrix<double> & K)
 {
-  // REPLACE WITH YOUR CODE
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
+  igl::min_quad_with_fixed_precompute(L, b, Eigen::SparseMatrix<double>(), false, data);
+
+  Eigen::MatrixXd C(F.rows(), 3);
+  igl::cotmatrix_entries(V, F, C);
+
+  std::vector<Eigen::Triplet<double>> triplets;
+
+  K.resize(V.rows(), 3 * V.rows());
+  for (int f = 0; f < F.rows(); f++) {
+    Eigen::RowVector3i face = F.row(f);
+    for (int x = 0; x < 3; x++) {
+      int i = face((x + 1) % 3);
+      int j = face((x + 2) % 3);
+      Eigen::RowVector3d edge = V.row(i) - V.row(j);
+      Eigen::Vector3d eij = C(f, x) / 3.0 * edge ;
+
+      for (int a = 0; a < 3; a++) {
+        int k = face((x + a) % 3);
+        for (int b = 0; b < 3; b++) {
+          triplets.emplace_back(Eigen::Triplet<double>(i, 3 * k + b, eij(b)));
+          triplets.emplace_back(Eigen::Triplet<double>(j, 3 * k + b, -eij(b)));
+        }
+      }
+    }
+  }
+  K.setFromTriplets(triplets.begin(), triplets.end());
+
 }

src/arap_single_iteration.cpp

@@ -1,10 +1,25 @@
+#include <igl/min_quad_with_fixed.h>
 #include "arap_single_iteration.h"
+#include <igl/polar_svd3x3.h>
 
 void arap_single_iteration(
   const igl::min_quad_with_fixed_data<double> & data,
   const Eigen::SparseMatrix<double> & K,
   const Eigen::MatrixXd & bc,
   Eigen::MatrixXd & U)
 {
-  // REPLACE WITH YOUR CODE
+  Eigen::MatrixXd C = K.transpose() * U;
+  Eigen::MatrixXd R(3 * U.rows(), 3);
+
+  for (int k = 0; k < U.rows(); k++) {
+    Eigen::Matrix3d Ck, Rk;
+    for (int i = 0; i < 3; i++) {
+      Ck.row(i) = C.row(3 * k + i);
+    }
+    igl::polar_svd3x3(Ck, Rk);
+    for (int i = 0; i < 3; i++) {
+      R.row(3 * k + i) = Rk.row(i);
+    }
+  }
+  igl::min_quad_with_fixed_solve(data, K * R, bc, Eigen::MatrixXd(), U);
 }

src/biharmonic_precompute.cpp

@@ -1,13 +1,24 @@
 #include "biharmonic_precompute.h"
 #include <igl/min_quad_with_fixed.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
 
 void biharmonic_precompute(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   const Eigen::VectorXi & b,
   igl::min_quad_with_fixed_data<double> & data)
 {
-  // REPLACE WITH YOUR CODE
-  data.n = V.rows();
+  Eigen::SparseMatrix<double> L, M;
+  igl::cotmatrix(V, F, L);
+  igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_DEFAULT, M);
+
+  Eigen::SparseMatrix<double> M_i(M.rows(), M.cols());
+  for (int i = 0; i < M.rows(); i++) {
+    M_i.coeffRef(i, i) = 1.0 / M.coeff(i, i);
+  }
+
+  Eigen::SparseMatrix<double> Q = L.transpose() * M_i * L;
+  igl::min_quad_with_fixed_precompute(Q, b, Eigen::SparseMatrix<double>(), false, data);
 }
 

src/biharmonic_solve.cpp

@@ -6,7 +6,7 @@ void biharmonic_solve(
   const Eigen::MatrixXd & bc,
   Eigen::MatrixXd & D)
 {
-  // REPLACE WITH YOUR CODE
-  D = Eigen::MatrixXd::Zero(data.n,3);
+  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(data.n, 1);
+  igl::min_quad_with_fixed_solve(data, B, bc, Eigen::MatrixXd(), D);
 }