From 8e4934fdc462c5b13e805cae3942ebfc6a3caab9 Mon Sep 17 00:00:00 2001
From: Zhiyuan <jerryz.cheng@mail.utoronto.ca>
Date: Mon, 15 Oct 2018 14:34:19 -0400
Subject: [PATCH] hw03 submit

---
 src/cotmatrix.cpp  | 42 ++++++++++++++++++++++++++++++++++++++----
 src/massmatrix.cpp | 17 +++++++++++++----
 src/smooth.cpp     | 39 ++++++++++++++++++++++++++++++++-------
 3 files changed, 83 insertions(+), 15 deletions(-)

diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..06b2a57 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,10 +1,44 @@
 #include "cotmatrix.h"
+#include <math.h>
+
+double cot(double a, double b, double c) {
+	double s = (a + b + c) / 2.0;
+	double A = sqrt(s * (s - a) * (s - b) * (s - c));
+	return (b*b + c*c - a*a) / (4.0 * A);
+}
+
 
 void cotmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::SparseMatrix<double> & L)
+	const Eigen::MatrixXd & l,
+	const Eigen::MatrixXi & F,
+	Eigen::SparseMatrix<double> & L)
 {
-  // Add your code here
+	int n = F.maxCoeff()+1;
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	tripletList.reserve(9 * F.rows());
+	
+	for (int i = 0; i < F.rows(); i++) {
+		double a = l(i, 0);
+		double b = l(i, 1);
+		double c = l(i, 2);
+		// case i == j
+		tripletList.push_back(T(F(i, 0), F(i, 0), -cot(b, c, a)/2.0));
+		tripletList.push_back(T(F(i, 0), F(i, 0), -cot(c, a, b)/2.0));
+		tripletList.push_back(T(F(i, 1), F(i, 1), -cot(a, b, c)/2.0));
+		tripletList.push_back(T(F(i, 1), F(i, 1), -cot(c, a, b)/2.0));
+		tripletList.push_back(T(F(i, 2), F(i, 2), -cot(a, b, c)/2.0));
+		tripletList.push_back(T(F(i, 2), F(i, 2), -cot(b, c, a)/2.0));
+		// case i != j
+		tripletList.push_back(T(F(i, 0), F(i, 2), cot(b, c, a)/2.0));
+		tripletList.push_back(T(F(i, 0), F(i, 1), cot(c, a, b)/2.0));
+		tripletList.push_back(T(F(i, 1), F(i, 2), cot(a, b, c)/2.0));
+		tripletList.push_back(T(F(i, 1), F(i, 0), cot(c, a, b)/2.0));
+		tripletList.push_back(T(F(i, 2), F(i, 1), cot(a, b, c)/2.0));
+		tripletList.push_back(T(F(i, 2), F(i, 0), cot(b, c, a)/2.0));
+	}
+
+	L.resize(n, n);
+	L.setFromTriplets(tripletList.begin(), tripletList.end());
 }
 
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..f10bda3 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,19 @@
 #include "massmatrix.h"
+#include <igl/doublearea.h>
 
 void massmatrix(
-  const Eigen::MatrixXd & l,
-  const Eigen::MatrixXi & F,
-  Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
+	const Eigen::MatrixXd & l,
+	const Eigen::MatrixXi & F,
+	Eigen::DiagonalMatrix<double, Eigen::Dynamic> & M)
 {
-  // Add your code here
+	int n = F.maxCoeff()+1;
+	M.resize(n);
+	Eigen::MatrixXd A;
+	igl::doublearea(l, A);
+	for (int i = 0; i < F.rows(); i++) {
+		for (int j = 0; j < 3; j++) {
+			M.diagonal()[F(i, j)] += A(i) / 6.0;
+		}
+	}
 }
 
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..f0a45ae 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,12 +1,37 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <igl/edge_lengths.h>
+#include <Eigen/SparseCholesky>
 
 void smooth(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXd & G,
-    double lambda,
-    Eigen::MatrixXd & U)
+	const Eigen::MatrixXd & V,
+	const Eigen::MatrixXi & F,
+	const Eigen::MatrixXd & G,
+	double lambda,
+	Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = G;
+	// get L and M
+	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);
+
+	// convert M from DiagonalMatrix to SparseMatrix
+		/* Eigen::SparseMatrix<double>(M) not work. */
+	typedef Eigen::Triplet<double> T;
+	std::vector<T> tripletList;
+	tripletList.reserve(M.rows());
+
+	for (int i = 0; i < M.rows(); i++) {tripletList.push_back(T(i, i, M.diagonal()[i]));}
+	Eigen::SparseMatrix<double> MS;
+	MS.resize(M.rows(), M.rows());
+	MS.setFromTriplets(tripletList.begin(), tripletList.end());
+
+	// compute U_{t+1}
+	U.resizeLike(G);
+	Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> cholesky(MS - lambda * L);
+	U = cholesky.solve(M * G);
 }