From 7559747e15967ea2e51570a474b4419a0a0066a8 Mon Sep 17 00:00:00 2001
From: eriszhang <jiayieris.zhang@mail.utoronto.ca>
Date: Sat, 13 Oct 2018 17:18:51 -0400
Subject: [PATCH 1/2] it works

---
 src/cotmatrix.cpp  | 52 ++++++++++++++++++++++++++++++++++++++++++++--
 src/massmatrix.cpp | 15 +++++++++++--
 src/smooth.cpp     | 25 ++++++++++++++++++++--
 3 files changed, 86 insertions(+), 6 deletions(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..cbeb4cd 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,58 @@
 #include "cotmatrix.h"
+#include <cmath>
+
+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;
+  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); // AC
+    double c = l(i, 2); // CA
+    // 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(F.maxCoeff() + 1, F.maxCoeff() + 1);
+  L.setFromTriplets(triplets.begin(), triplets.end());
 }
-
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..a9bd566 100644
--- a/src/massmatrix.cpp
+++ b/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;
+    }
+  }
 }
-
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..30a2618 100644
--- a/src/smooth.cpp
+++ b/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);
 }

From fad4683a12f6c6ab098d102770989ac1085ce57c Mon Sep 17 00:00:00 2001
From: eriszhang <jiayieris.zhang@mail.utoronto.ca>
Date: Sun, 14 Oct 2018 22:39:17 -0400
Subject: [PATCH 2/2] fix a typo

---
 src/cotmatrix.cpp | 12 +++++++++---
 src/smooth.cpp    |  2 +-
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index cbeb4cd..93cfb53 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,5 +1,6 @@
 #include "cotmatrix.h"
 #include <cmath>
+#include <iostream>
 
 using namespace Eigen;
 using namespace std;
@@ -25,6 +26,10 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   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);
@@ -32,8 +37,8 @@ void cotmatrix(
     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); // AC
-    double c = l(i, 2); // CA
+    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);
@@ -53,6 +58,7 @@ void cotmatrix(
     triplets.push_back(Triplet<double>(VA, VA, -cotB / 2.0));
     triplets.push_back(Triplet<double>(VC, VC, -cotB / 2.0));
   }
-  L.resize(F.maxCoeff() + 1, F.maxCoeff() + 1);
+
+  L.resize(n, n);
   L.setFromTriplets(triplets.begin(), triplets.end());
 }
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 30a2618..0e27eb8 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -30,4 +30,4 @@ void smooth(
   SimplicialLDLT<SparseMatrix<double>> solver(right);
   MatrixXd left = M * G;
   U = solver.solve(left);
-}
+}
\ No newline at end of file