From 0f1372fa48d9cc37282f1666029146195aff9ee0 Mon Sep 17 00:00:00 2001
From: augustin64 <me.git@augustin64.fr>
Date: Wed, 19 Mar 2025 12:15:26 +0100
Subject: [PATCH] Try to do Q6

---
 .../c++/dirichletProblems.cpp                 | 75 ++++++++++++++++++-
 1 file changed, 73 insertions(+), 2 deletions(-)

diff --git a/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp b/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp
index 670fbec..d22b9a3 100644
--- a/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp
+++ b/geometry-processing/6-GeomProcessing/c++/dirichletProblems.cpp
@@ -77,21 +77,24 @@ int main(int argc, char **argv) {
                                  geometry->inputVertexPositions,
                                  mesh->getFaceVertexList(),
                                  polyscopePermutations(*mesh));
+  geometry->requireEdgeCotanWeights();
   geometry->requireVertexPositions();
   geometry->requireFaceNormals();
   geometry->requireFaceAreas();
 
+  auto vertCount = geometry->inputVertexPositions.size();
+
   std::vector<double> sin_coords;
   {//* First question : showing basic quantity
     std::vector<double> snd_coordinate;
-    for (auto i=0; i < geometry->inputVertexPositions.size(); i++) {
+    for (auto i=0; i < vertCount; i++) {
       auto pos = geometry->inputVertexPositions[i];
       snd_coordinate.push_back((double)pos[1]);
     }
     sfm->addVertexScalarQuantity("Second coordinate", snd_coordinate);
 
     float lambda = 18.88;
-    for (auto i=0; i < geometry->inputVertexPositions.size(); i++) {
+    for (auto i=0; i < vertCount; i++) {
       auto pos = geometry->inputVertexPositions[i];
       double posX = cos(PI / 4) * pos.x - sin(PI / 4) * pos.y;
       double posY = sin(PI / 4) * pos.x + cos(PI / 4) * pos.y;
@@ -124,6 +127,74 @@ int main(int argc, char **argv) {
   }
 
   //* Laplace-Beltrami
+  {//* 6th question : implement the Laplace-Beltrami operator on triangular meshes
+    //* Create Laplace-Beltrami in a standard matrix
+    Eigen::SparseMatrix<double> laplaceBeltrami(vertCount, vertCount);
+
+    {//* Create LaplaceBeltrami Lij
+      // We create the Lii and Lij (i!=j) values in different ways : Lii are added to the triplet list at the end
+      std::vector<Eigen::Triplet<double>> LtripletList;
+      std::vector<double> Lii (vertCount, 0); // for Lii values
+
+      for (Vertex v1 : mesh->vertices()) {
+        for (Edge e : v1.adjacentEdges()) {
+          auto wij = geometry->edgeCotanWeight(e);
+
+          Lii.at(v1.getIndex()) += wij;
+
+          auto v2 = e.firstVertex().getIndex() == v1.getIndex() ? e.secondVertex() : e.firstVertex();
+          LtripletList.push_back(Eigen::Triplet<double>(v1.getIndex(), v2.getIndex(), -wij));
+
+          assert(v1.getIndex() != v2.getIndex());
+        }
+      }
+      
+      for (Vertex v : mesh->vertices()) { // add Lii triplets
+        LtripletList.push_back(Eigen::Triplet<double>(v.getIndex(), v.getIndex(), Lii.at(v.getIndex())));
+      }
+      laplaceBeltrami.setFromTriplets(LtripletList.begin(), LtripletList.end());
+      laplaceBeltrami.finalize();
+
+      //TODO: check that row and columns sum to 0
+      //double sumCol, sumRow;
+      //for (Vertex v1: mesh->vertices()) {
+      //  sumCol = sumRow = 0;
+      //  for (Vertex v2: mesh->vertices()) {
+      //    sumCol += laplaceBeltrami((size_t)v1.getIndex(), (size_t)v2.getIndex()); // does not work
+      //    sumRow += laplaceBeltrami((size_t)v2.getIndex(), (size_t)v1.getIndex());
+      //  }
+      //  std::cout << sumCol << " " << sumRow << std::endl;
+      //}
+    }
+
+    Eigen::SparseMatrix<double> M_inv(vertCount, vertCount);
+    {//* Create M in a sparse matrix
+      std::vector<Eigen::Triplet<double>> MtripletList;
+
+      for (Vertex v : mesh->vertices()) {
+        double val = 0;
+        for (auto f : v.adjacentFaces()) {
+          val += 1/3*geometry->faceArea(f);
+        }
+        assert(val != 0);
+        MtripletList.push_back(Eigen::Triplet<double>(v.getIndex(), v.getIndex(), 1/val));
+      }
+      M_inv.setFromTriplets(MtripletList.begin(), MtripletList.end());
+      M_inv.finalize();
+    }
+    
+    //* Show the result on some function of q1
+    Vector<double> eigen_sin_coords(vertCount); // We need to convert this to an Eigen::Vector
+    for (auto i=0; i < vertCount; i++) {
+      eigen_sin_coords[i] = sin_coords.at(i);
+    }
+
+    Eigen::SparseMatrix<double> DeltaCotan = M_inv * laplaceBeltrami;
+    geometrycentral::Vector<double> lb_sin_coords = DeltaCotan * eigen_sin_coords;
+
+    sfm->addVertexScalarQuantity("Laplace-Beltrami of Sinusoid data", lb_sin_coords);
+  }
+
 
   // Give control to the polyscope gui
   polyscope::show();