From 5e3acf4287638bdc8e05d2b2c956050e5f88d8cd Mon Sep 17 00:00:00 2001
From: Simeon Krastnikov <simeon.krastnikov@mail.utoronto.ca>
Date: Thu, 11 Oct 2018 22:14:47 -0400
Subject: [PATCH] HW3 Submission

---
 src/cotmatrix.cpp  | 27 +++++++++++++++++++++++++--
 src/massmatrix.cpp | 12 ++++++++++--
 src/smooth.cpp     | 28 ++++++++++++++++++++++++++--
 3 files changed, 61 insertions(+), 6 deletions(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..9524d5f 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,33 @@
 #include "cotmatrix.h"
+#include <vector>
+#include <math.h>
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tripletList;
+  
+  for (int f = 0; f < F.rows(); f++) {
+    for (int v = 0; v < 3; v++) {
+      int i = F(f, v);
+      int j = F(f, (v+1) % 3);
+      
+      double eij = l(f, v);
+      double ej = l(f, (v+1) % 3);
+      double ei = l(f, (v+2) % 3);
+      
+      double s = (ei + ej + eij) / 2;
+      double A = sqrt(s * (s-ei) * (s-ej) * (s-eij));  // heron's formula
+      double cot_aij = (ei*ei + ej*ej - eij*eij) / (4 * A);  // law of cosines + sine area formula
+      tripletList.push_back(T(i, j, cot_aij / 2));
+      tripletList.push_back(T(j, i, cot_aij / 2));
+      tripletList.push_back(T(i, i, -cot_aij / 2));
+      tripletList.push_back(T(j, j, -cot_aij / 2));
+    }
+  }
+  
+  L.setFromTriplets(tripletList.begin(), tripletList.end());
 }
-
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..1cb6ec5 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,18 @@
 #include "massmatrix.h"
+#include <math.h>
 
 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);
+  for (int f = 0; f < F.rows(); f++) {
+    double s = l.row(f).sum() / 2;
+    double f_A = sqrt(s * (s-l(f,0)) * (s-l(f,1)) * (s-l(f,2)));
+    for (int v = 0; v < 3; v++) {
+      int i = F(f, v);
+      M.diagonal()[i] += f_A / 3;
+    }
+  }
 }
-
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..e8e7145 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,9 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include "igl/edge_lengths.h"
+#include <vector>
+#include <Eigen/SparseCholesky>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -7,6 +12,25 @@ void smooth(
     double lambda,
     Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+  Eigen::MatrixXd l(F.rows(), 3);
+  igl::edge_lengths(V, F, l);
+
+  int n = F.maxCoeff() + 1;
+  Eigen::SparseMatrix<double> L(n,n);
+  Eigen::DiagonalMatrix<double,Eigen::Dynamic> M;
+  cotmatrix(l, F, L);
+  massmatrix(l, F, M);
+ 
+  // construct the system matrix A = M - lambda * L
+  Eigen::SparseMatrix<double> A(n,n);
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tripletList;
+  for (int i = 0; i < M.rows(); i++) {
+    tripletList.push_back(T(i, i, M.diagonal()[i]));
+  }
+  A.setFromTriplets(tripletList.begin(), tripletList.end());
+  A -= lambda * L;
+  
+  Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
+  U = chol.solve(M*G);
 }