diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..2d5519f 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,4 +1,11 @@
 #include "lscm.h"
+#include "vector_area_matrix.h"
+#include <Eigen/SVD>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/eigs.h>
+#include <igl/repdiag.h>
 
 void lscm(
   const Eigen::MatrixXd & V,
@@ -6,5 +13,25 @@ void lscm(
   Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = V.leftCols(2);
+  // U = V.leftCols(2);
+  Eigen::SparseMatrix<double> L,A,Q;
+  igl::cotmatrix(V,F,L);
+  vector_area_matrix(F,A);
+  igl::repdiag(L,2,Q);
+  Q = Q - A;
+  Eigen::SparseMatrix<double> M,B;
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
+  igl::repdiag(M,2,B);
+
+  Eigen::MatrixXd sU;
+  Eigen::VectorXd sS;
+  igl::eigs(Q,B,3,igl::EIGS_TYPE_SM,sU,sS);
+
+  U.resize(V.rows(), 2);
+  U.col(0) = sU.col(2).head(V.rows());
+  U.col(1) = sU.col(2).tail(V.rows());
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(U.transpose() * U, Eigen::ComputeThinU | Eigen::ComputeThinV);
+  U = U * svd.matrixU();
+  // since U is flipped due to boundary issue, I times -1 to flip the normal.
+  U.col(0) *= -1;
 }
diff --git a/src/tutte.cpp b/src/tutte.cpp
index 6f6109b..9215a58 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,4 +1,8 @@
 #include "tutte.h"
+#include <igl/min_quad_with_fixed.h>
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/edges.h>
 
 void tutte(
   const Eigen::MatrixXd & V,
@@ -6,6 +10,33 @@ void tutte(
   Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = V.leftCols(2);
-}
+  // U = V.leftCols(2);
+  Eigen::VectorXi bnd;
+  igl::boundary_loop(F,bnd);
 
+  Eigen::MatrixXi E;
+  igl::edges(F,E);
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tlist;
+  tlist.reserve(E.rows()*2);
+  for (int i=0; i<E.rows();i++) {
+    tlist.push_back(T(E(i,0),E(i,1),1.0));
+    tlist.push_back(T(E(i,1),E(i,0),1.0));
+  }
+  Eigen::SparseMatrix<double> L;
+  int L_size = E.maxCoeff() + 1;
+  L.resize(L_size,L_size);
+  L.setFromTriplets(tlist.begin(),tlist.end());
+
+  for(int i=0; i<L_size; i++) {
+    L.insert(i,i) = -L.row(i).sum(); 
+  }
+
+  Eigen::MatrixXd UV;
+  igl::map_vertices_to_circle(V,bnd,UV);
+
+  Eigen::MatrixXd B = Eigen::MatrixXd::Zero(V.rows(),2);
+  igl::min_quad_with_fixed_data<double> data;
+  igl::min_quad_with_fixed_precompute(L,bnd,Eigen::SparseMatrix<double>(),false,data);
+  igl::min_quad_with_fixed_solve(data,B,UV,Eigen::MatrixXd(),U);
+}
\ No newline at end of file
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..ab49b31 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,11 +1,38 @@
 #include "vector_area_matrix.h"
+#include <igl/boundary_loop.h>
 
 void vector_area_matrix(
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double>& A)
 {
   // Replace with your code
-  int V_size = F.maxCoeff()+1;
-  A.resize(V_size*2,V_size*2);
+   int V_size = F.maxCoeff()+1;
+   A.resize(V_size*2,V_size*2);
+
+  std::vector<std::vector<int>> bnd;
+  igl::boundary_loop(F,bnd);
+
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> tlist;
+
+  for (int i=0; i<bnd.size(); i++) { // row
+    for (int j=0; j<bnd[i].size(); j++) { // col
+      if (j == bnd[i].size() - 1){
+        tlist.push_back(T(bnd[i][j],bnd[i][0] + V_size,0.5));
+        tlist.push_back(T(bnd[i][0], bnd[i][j] + V_size,-0.5));
+        tlist.push_back(T(bnd[i][0] + V_size,bnd[i][j], 0.5));
+        tlist.push_back(T(bnd[i][j] + V_size,bnd[i][0], -0.5));
+      }
+      else {
+        tlist.push_back(T(bnd[i][j],bnd[i][j+1] + V_size,0.5));
+        tlist.push_back(T(bnd[i][j+1], bnd[i][j] + V_size,-0.5));
+        tlist.push_back(T(bnd[i][j+1]+V_size,bnd[i][j], 0.5));
+        tlist.push_back(T(bnd[i][j] + V_size,bnd[i][j+1], -0.5));
+      }
+    }
+  }
+  A.setFromTriplets(tlist.begin(),tlist.end());
+
+
 }