From 4332d115aa562490c6a6ad4256ddc9af5074fd30 Mon Sep 17 00:00:00 2001
From: Meng-Wei <mmeng.wei@mail.utoronto.ca>
Date: Mon, 15 Oct 2018 17:23:32 -0400
Subject: [PATCH] final version

---
 src/cotmatrix.cpp  | 29 ++++++++++++++++++++++++++++-
 src/massmatrix.cpp | 18 +++++++++++++++++-
 src/smooth.cpp     | 27 +++++++++++++++++++++++++--
 3 files changed, 70 insertions(+), 4 deletions(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..74d2dee 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,37 @@
 #include "cotmatrix.h"
+#include <iostream>
+
+typedef Eigen::Triplet<double> T;
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+  // setFromTriplets will sum up all entries:
+  std::vector<T> triplet_list;
+  for (int n = 0; n < F.rows(); n ++) {
+    for (int m = 0; m < 3; m ++) {
+      int i = F(n, m);
+      int j = F(n, (m+1)%3);
+      int k = F(n, (m+2)%3);
+
+      double ei = l(n, m);
+      double ej = l(n, (m+1)%3);
+      double ek = l(n, (m+2)%3);
+
+      double cos_ij = (pow(ei, 2) + pow(ej, 2) - pow(ek, 2)) / (2.0*ei*ej);
+      double sin_ij = sqrt(1-pow(cos_ij, 2));
+      double cot_ij = 0.5 * (cos_ij / sin_ij);
+
+      triplet_list.push_back(T(i, j, cot_ij));
+      triplet_list.push_back(T(j, i, cot_ij));
+      triplet_list.push_back(T(i, i, -cot_ij));
+      triplet_list.push_back(T(j, j, -cot_ij));
+    }
+  }
+  int num_V = F.maxCoeff() + 1;
+  L.resize(num_V, num_V);
+  L.setFromTriplets(triplet_list.begin(), triplet_list.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..85401d2 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,26 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
+#include <iostream>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
+  // Calculate Area
+  Eigen::MatrixXd dblA;
+  igl::doublearea(l, dblA);
+  dblA = dblA / 2.0;
+
+  // Constuct M:
+  int num_V = F.maxCoeff() + 1;
+  M.setZero(num_V);
+  for(int i = 0; i < F.rows(); i ++){
+    double area = dblA(i);
+    for(int j = 0; j < 3; j ++) {
+      int cur = F(i, j);
+      M.diagonal()(cur) += area / 3.0;
+    }
+  }
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..8d219af 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,5 +1,10 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <igl/edge_lengths.h>
+#include <iostream>
 
+typedef Eigen::Triplet<double> T;
 void smooth(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
@@ -7,6 +12,24 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+  // Find Laplacian
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+  Eigen::SparseMatrix<double> L;
+  cotmatrix(l, F, L);
+
+  // Find M
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+  massmatrix(l, F, M);
+  std::vector<T> triplet_list;
+  for(int i = 0; i < M.rows(); i ++) {
+    triplet_list.push_back(T(i, i, M.diagonal()(i)));
+  }
+  Eigen::SparseMatrix<double> sparse_M(M.rows(), M.rows());
+  sparse_M.setFromTriplets(triplet_list.begin(), triplet_list.end());
+
+  // Solve system:
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+  solver.compute(-lambda * L + sparse_M);
+  U = solver.solve(M * G);
 }