Jiayi(Eris) Zhang HW3

src/cotmatrix.cpp

@@ -1,10 +1,64 @@
 #include "cotmatrix.h"
+#include <cmath>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+double compute_cot(
+  double a,
+  double b,
+  double c)
+{
+  // compute area: Heron's formula
+  double s = (a + b + c) / 2.0;
+  double area = sqrt(s * (s - a) * (s - b) * (s - c));
+  // cotA = cosA / sinA
+  double sinA = 2.0 * area / (b * c);
+  double cosA = (b * b + c * c - a * a) * 1.0 / (2 * b * c);
+  double cotA = cosA / sinA;
+  return cotA;
+}
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
-}
+  vector<Triplet<double>> triplets;
+  int n = F.maxCoeff() + 1;
+  triplets.reserve(n * 12);
+
+
+  for (int i = 0; i < F.rows(); i++) {
+    // vertices
+    int VA = F(i, 0);
+    int VB = F(i, 1);
+    int VC = F(i, 2);
+    // edges l: columns correspond to edges 23,31,12
+    double a = l(i, 0); // BC
+    double b = l(i, 1); // CA
+    double c = l(i, 2); // AB
+    // cot
+    double cotA = compute_cot(a, b, c);
+    double cotB = compute_cot(b, c, a);
+    double cotC = compute_cot(c, a, b);
+    // symmetric
+    triplets.push_back(Triplet<double>(VA, VB, cotC / 2.0));
+    triplets.push_back(Triplet<double>(VB, VA, cotC / 2.0));
+    triplets.push_back(Triplet<double>(VB, VC, cotA / 2.0));
+    triplets.push_back(Triplet<double>(VC, VB, cotA / 2.0));
+    triplets.push_back(Triplet<double>(VC, VA, cotB / 2.0));
+    triplets.push_back(Triplet<double>(VA, VC, cotB / 2.0));
+    // diagonal
+    triplets.push_back(Triplet<double>(VA, VA, -cotC / 2.0));
+    triplets.push_back(Triplet<double>(VB, VB, -cotC / 2.0));
+    triplets.push_back(Triplet<double>(VB, VB, -cotA / 2.0));
+    triplets.push_back(Triplet<double>(VC, VC, -cotA / 2.0));
+    triplets.push_back(Triplet<double>(VA, VA, -cotB / 2.0));
+    triplets.push_back(Triplet<double>(VC, VC, -cotB / 2.0));
+  }
 
+  L.resize(n, n);
+  L.setFromTriplets(triplets.begin(), triplets.end());
+}

src/massmatrix.cpp

@@ -1,10 +1,21 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
+
+using namespace Eigen;
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
+  M.resize(F.maxCoeff() + 1);
+  MatrixXd doubleArea;
+  igl::doublearea(l, doubleArea);
+  for (int i = 0; i < F.rows(); i++) {
+    for (int j = 0; j < 3; j++) {
+      int v = F(i, j);
+      double area = 0.5 * doubleArea(i, 0);
+      M.diagonal()[v] += area / 3.0;
+    }
+  }
 }
-

src/smooth.cpp

@@ -1,4 +1,10 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <Eigen/SparseCholesky>
+#include <igl/edge_lengths.h>
+
+using namespace Eigen;
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -7,6 +13,21 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
-}
+  // construct cotmatrix and massMatrix
+  MatrixXd l;
+  igl::edge_lengths(V, F, l);
+  SparseMatrix<double> L;
+  DiagonalMatrix<double, Dynamic> M;
+  cotmatrix(l, F, L);
+  massmatrix(l, F, M);
+
+  // solve linear equation: M * G = (M - lam * L) * U
+  SparseMatrix<double> right(M.rows(), M.rows());
+  for(int i = 0; i < M.rows(); i++) {
+  	right.coeffRef(i, i) += M.diagonal()[i];
+  }
+  right -= lambda * L;
+  SimplicialLDLT<SparseMatrix<double>> solver(right);
+  MatrixXd left = M * G;
+  U = solver.solve(left);
+}