Construct cotan matrix and mass matrix

Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h

@@ -30,18 +30,25 @@
 #include <CGAL/property_map.h>
 
 #include <Eigen/Cholesky>
-
+#include <Eigen/Sparse>
-#include <Eigen/Dense>
 
 #include <boost/foreach.hpp>
-
+#include <CGAL/boost/graph/properties.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Dynamic_property_map.h>
+#include <vector>
+#include <CGAL/Vector_3.h>
+#include <CGAL/squared_distance_3.h>
+#include <CGAL/Polygon_mesh_processing/measure.h>
 namespace CGAL {
 namespace Heat_method_3 {
 
   // This class will later go into another file
   // It encapsulates what we use from Eigen so that one potentially can use another LA library
   struct Heat_method_Eigen_traits_3 {
-    typedef Eigen::Matrix3d Matrix;
+    typedef Eigen::SparseMatrix<double> SparseMatrix;
+    typedef Eigen::Triplet<double> T;
   };
 
 
@@ -67,11 +74,13 @@ namespace Heat_method_3 {
     typedef typename graph_traits::edge_descriptor                edge_descriptor;
     typedef typename graph_traits::halfedge_descriptor        halfedge_descriptor;
     typedef typename graph_traits::face_descriptor                face_descriptor;
-
+    typedef typename graph_traits::vertex_iterator                vertex_iterator;
     /// Geometric typedefs
     typedef typename Traits::Point_3                                      Point_3;
     typedef typename Traits::FT                                                FT;
 
+    typedef typename Simple_cartesian<double>::Vector_3                    vector;
+
     // The Vertex_id_map is a property map where you can associate an index to a vertex_descriptor
     typedef typename boost::graph_traits<TriangleMesh>::vertices_size_type vertices_size_type;
 
@@ -79,7 +88,8 @@ namespace Heat_method_3 {
     typedef typename boost::property_map<TriangleMesh, Vertex_property_tag >::type Vertex_id_map;
     Vertex_id_map vertex_id_map;
 
-    typedef typename LA::Matrix Matrix;
+    typedef typename LA::SparseMatrix Matrix;
+    typedef typename LA::T triplet;
   public:
 
     Heat_method_3(const TriangleMesh& tm_)
@@ -121,7 +131,7 @@ namespace Heat_method_3 {
     /**
      *return current source set
     */
-    vertex_descriptor getSources()
+    vertex_descriptor get_sources()
     {
       return sources;
     }
@@ -129,7 +139,7 @@ namespace Heat_method_3 {
     /**
      * clear the current source set
      */
-    void clearSources()
+    void clear_sources()
     {
       sources.clear();
       return;
@@ -138,14 +148,14 @@ namespace Heat_method_3 {
     /**
      * return vertex_descriptor to first vertex in the source set
      */
-    vertex_descriptor sourcesBegin()
+    vertex_iterator sources_begin()
     {
       return sources.begin();
     }
     /**
      * return vertex_descriptor to last vertex in the source set
      */
-    vertex_descriptor sourcesEnd()
+    vertex_descriptor sources_end()
     {
       return sources.end();
     }
@@ -173,6 +183,64 @@ namespace Heat_method_3 {
           std::cout << p << std::endl;
         }
       }
+
+      int m = num_vertices(tm);
+      Matrix c(m,m);
+      Matrix A(m,m);
+      std::vector<triplet> A_matrix_entries;
+      std::vector<triplet> c_matrix_entries;
+      CGAL::Vertex_around_face_iterator<TriangleMesh> vbegin, vend;
+      BOOST_FOREACH(face_descriptor f, faces(tm)) {
+        boost::tie(vbegin, vend) = vertices_around_face(halfedge(f,tm),tm);
+        vertex_descriptor current = *(vbegin);
+        vertex_descriptor neighbor_one = *(vbegin++);
+        vertex_descriptor neighbor_two = *(vend);
+        double i = get(vertex_id_map, current);
+        double j = get(vertex_id_map, neighbor_one);
+        double k = get(vertex_id_map, neighbor_two);
+        Point_3 p_i = get(vpm,current);
+        Point_3 p_j = get(vpm, neighbor_one);
+        Point_3 p_k = get(vpm, neighbor_two);
+        //If the passed in mesh is not a triangle mesh, the algorithm breaks here
+        vector cross = CGAL::cross_product((p_j-p_i), (p_k-p_i));
+        double dot = to_double((p_j-p_i)*(p_k-p_i));
+        double squared_cross = to_double(CGAL::approximate_sqrt(cross*cross));
+        double cotan_i = dot/squared_cross;
+        c_matrix_entries.push_back(triplet(j,k ,-.5*cotan_i));
+        c_matrix_entries.push_back(triplet(k,j,-.5* cotan_i));
+        c_matrix_entries.push_back(triplet(j,j,.5* cotan_i));
+        c_matrix_entries.push_back(triplet(k,k,.5* cotan_i));
+
+        cross = CGAL::cross_product((p_i-p_j), (p_k-p_j));
+        dot = to_double((p_i-p_j)*(p_k-p_j));
+        squared_cross = to_double(CGAL::approximate_sqrt(cross*cross));
+        double cotan_j = dot/squared_cross;
+        c_matrix_entries.push_back(triplet(i,k ,-.5*cotan_j));
+        c_matrix_entries.push_back(triplet(k,i,-.5* cotan_j));
+        c_matrix_entries.push_back(triplet(i,i,.5* cotan_j));
+        c_matrix_entries.push_back(triplet(k,k,.5* cotan_j));
+
+        cross = CGAL::cross_product((p_i-p_k), (p_j-p_k));
+        dot = to_double((p_i-p_k)*(p_j-p_k));
+        squared_cross = to_double(CGAL::approximate_sqrt(cross*cross));
+        double cotan_k = dot/squared_cross;
+        c_matrix_entries.push_back(triplet(i,j,-.5*cotan_k));
+        c_matrix_entries.push_back(triplet(j,i,-.5* cotan_k));
+        c_matrix_entries.push_back(triplet(i,i,.5* cotan_k));
+        c_matrix_entries.push_back(triplet(j,j,.5* cotan_k));
+
+        double area_face = CGAL::Polygon_mesh_processing::face_area(f,tm);
+        A_matrix_entries.push_back(triplet(i,i, (1/3)*area_face));
+        A_matrix_entries.push_back(triplet(j,j, (1/3)*area_face));
+        A_matrix_entries.push_back(triplet(k,k, (1/3)*area_face));
+      }
+
+      A.setFromTriplets(A_matrix_entries.begin(), A_matrix_entries.end());
+      c.setFromTriplets(c_matrix_entries.begin(), c_matrix_entries.end());
+
+
+
+
     }
 
     const TriangleMesh& tm;