Adrian She HW3 Submission

src/cotmatrix.cpp

@@ -1,10 +1,51 @@
 #include "cotmatrix.h"
+#include "igl/doublearea.h"
+#include <cmath>
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+
+  // Compute areas of each face
+  Eigen::VectorXd A;
+  igl::doublearea(l, A);
+
+  // Loop through faces to determine cotangent matrix elements
+  // Using Eigen's functionality that triplets at duplicate (i, j)
+  // are summed up when creating a sparse matrix
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> triples;
+  triples.reserve(F.rows() * 12);
+  for (int i = 0; i < F.rows(); i++) {
+        // Edge e = (F(j), F(j+1)) is opposite of F(j+2) mod 3
+	// cotangent of the opposite angle is (-e^2 + l_1^2 + l_2^2) / (4 A)
+	// where l_1, l_2 are the other sides and A is the area of current face. 
+	double len_01 = ((pow(l(i, 0), 2) + pow(l(i,1), 2) - pow(l(i,2), 2))) / (4.0 * A(i)); 
+	double len_12 = ((pow(l(i, 1), 2) + pow(l(i,2), 2) - pow(l(i,0), 2))) / (4.0 * A(i)); 
+	double len_02 = ((pow(l(i, 0), 2) + pow(l(i,2), 2) - pow(l(i,1), 2))) / (4.0 * A(i)); 
+        double entry_01 = len_01 * 0.5;
+        double entry_12 = len_12 * 0.5;
+  	double entry_02 = len_02 * 0.5;
+        // Insert entries
+        triples.push_back(T(F(i,0), F(i,1), entry_01)); 
+        triples.push_back(T(F(i,1), F(i,0), entry_01)); 
+        triples.push_back(T(F(i,0), F(i,0), -entry_01)); 
+        triples.push_back(T(F(i,1), F(i,1), -entry_01)); 
+
+        triples.push_back(T(F(i,0), F(i,2), entry_02)); 
+        triples.push_back(T(F(i,2), F(i,0), entry_02)); 
+        triples.push_back(T(F(i,0), F(i,0), -entry_02)); 
+        triples.push_back(T(F(i,2), F(i,2), -entry_02)); 
+
+        triples.push_back(T(F(i,1), F(i,2), entry_12)); 
+        triples.push_back(T(F(i,2), F(i,1), entry_12)); 
+        triples.push_back(T(F(i,1), F(i,1), -entry_12)); 
+        triples.push_back(T(F(i,2), F(i,2), -entry_12)); 
+  }
+  L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+  L.setFromTriplets(triples.begin(), triples.end()); 
 }
 

src/massmatrix.cpp

@@ -1,10 +1,27 @@
 #include "massmatrix.h"
+#include "igl/doublearea.h"
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
+
+  // Compute area of each face
+  Eigen::VectorXd A;
+  igl::doublearea(l, A);
+  M.resize(F.maxCoeff() + 1);
+
+  // Mass matrix entry (i, i) is 1/3 * (sum of faces containing vertex i)
+  for (int i = 0; i <= F.maxCoeff(); i++) {
+	double entry = 0;
+	for (int j = 0; j < F.rows(); j++) {
+		if (F(j, 0) == i || F(j,1) == i || F(j, 2) == i) {
+			entry = entry + A(j);
+		}
+	}
+	M.diagonal()[i] = entry / 3;
+  }
 }
 

src/smooth.cpp

@@ -1,4 +1,8 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include "igl/edge_lengths.h"
+#include "Eigen/SparseCholesky"
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +12,25 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+
+  // Compute edge lengths
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+
+  // Construct cot and mass matrices
+  Eigen::SparseMatrix<double> L;
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+  cotmatrix(l, F, L); 
+  massmatrix(l, F, M);
+
+  // Solve for U as solution of M G = (M - lambda * L) U;
+  Eigen::SparseMatrix<double> I(F.maxCoeff() + 1, F.maxCoeff() + 1);
+  I.setIdentity();
+  Eigen::SparseMatrix<double> L_lambda = I - lambda * M.inverse() * L;
+  L_lambda = M * L_lambda;
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+  solver.compute(L_lambda);
+  for (int i = 0; i < G.cols(); i++) {
+	U.col(i) = solver.solve(M * G.col(i));
+  }
 }