From 641684f50373e654b9f28af1721bb4cb90d5b486 Mon Sep 17 00:00:00 2001
From: Feilong Song <feilong.song@mail.utoronto.ca>
Date: Sun, 14 Oct 2018 16:44:37 -0400
Subject: [PATCH] works

---
 src/cotmatrix.cpp  | 35 +++++++++++++++++++++++++++++++++++
 src/massmatrix.cpp | 21 +++++++++++++++++++++
 src/smooth.cpp     | 24 +++++++++++++++++++++++-
 3 files changed, 79 insertions(+), 1 deletion(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..3772981 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -6,5 +6,40 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
+  int V_num = F.maxCoeff() + 1;
+  L.resize(V_num, V_num);
+
+  auto cot = [](double a, double b, double c){
+    double cosA = (b * b + c * c - a * a) * 1.0 / (2 * b * c);
+    double sinA = sqrt(1 - cosA * cosA);
+    return cosA / sinA;
+  };
+
+  for (int i = 0; i < F.rows(); i++) {
+    double edge0 = l(i, 0);
+    double edge1 = l(i, 1);
+    double edge2 = l(i, 2);
+
+    double cot0 = cot(edge0, edge1, edge2);
+    double cot1 = cot(edge1, edge0, edge2);
+    double cot2 = cot(edge2, edge0, edge1);
+
+    int v0 = F(i, 0), v1 = F(i, 1), v2 = F(i, 2);
+    L.coeffRef(v0, v1) += 0.5 * cot2;
+    L.coeffRef(v1, v0) += 0.5 * cot2;
+    L.coeffRef(v1, v2) += 0.5 * cot0;
+    L.coeffRef(v2, v1) += 0.5 * cot0;
+    L.coeffRef(v0, v2) += 0.5 * cot1;
+    L.coeffRef(v2, v0) += 0.5 * cot1;
+
+    L.coeffRef(v0, v0) -= 0.5 * cot1;
+    L.coeffRef(v0, v0) -= 0.5 * cot2;
+    L.coeffRef(v1, v1) -= 0.5 * cot0;
+    L.coeffRef(v1, v1) -= 0.5 * cot2;
+    L.coeffRef(v2, v2) -= 0.5 * cot0;
+    L.coeffRef(v2, v2) -= 0.5 * cot1;
+  }
+
+  L.makeCompressed();
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..1968c4a 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,4 +1,5 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +7,25 @@ void massmatrix(
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
   // Add your code here
+  int V_num = F.maxCoeff() + 1;
+  M.resize(V_num);
+
+  Eigen::VectorXd diag(V_num);
+
+  for (int i = 0; i < F.rows(); i++) {
+    double edge0 = l(i, 0);
+    double edge1 = l(i, 1);
+    double edge2 = l(i, 2);
+    double area = 0.25 * sqrt((edge0 + edge1 + edge2) * (-edge0 + edge1 + edge2) * (edge0 - edge1 + edge2) * (edge0 + edge1 - edge2));
+
+    int v0 = F(i, 0), v1 = F(i, 1), v2 = F(i, 2);
+
+    diag(v0) += area / 3;
+    diag(v1) += area / 3;
+    diag(v2) += area / 3;
+  }
+
+  M.diagonal() = diag;
+
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..2c66ace 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,8 @@
 #include "smooth.h"
+#include "massmatrix.h"
+#include "cotmatrix.h"
+#include <igl/edge_lengths.h>
+#include <Eigen/SparseCholesky>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +12,23 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V, F, l);
+
+  Eigen::SparseMatrix<double> L;
+  cotmatrix(l, F, L);
+
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M;
+  massmatrix(l, F, M);
+
+  Eigen::SparseMatrix<double> A = - lambda * L;
+  for (int i = 0; i < M.rows(); i++) {
+    A.coeffRef(i, i) += M.diagonal()(i);
+  }
+
+  Eigen::MatrixXd B = M * G;
+
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+  solver.compute(A);
+  U = solver.solve(B);
 }