From 40deb94656d8f34e9fc0aadb9baec6828ac4d7d0 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Tue, 17 Jul 2018 13:38:39 +0200
Subject: [PATCH] copied some fixes from Christinas branch

---
 .../include/CGAL/Heat_method_3/Heat_method_3.h   |  4 ++--
 .../Intrinsic_Delaunay_Triangulation_3.h         |  3 ++-
 .../heat_method_intrinsic_surface_mesh_test.cpp  | 16 +++++++++-------
 3 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
index 02d65bb7b7b..ab1b812960b 100644
--- a/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
+++ b/Heat_method_3/include/CGAL/Heat_method_3/Heat_method_3.h
@@ -260,7 +260,7 @@ namespace Heat_method_3 {
     {
       Eigen::VectorXd u;
       Matrix A = (M+ a_time_step*c);
-      Eigen::SimplicialLLT<Matrix> solver;
+      Eigen::SimplicialLDLT<Matrix> solver;
       solver.compute(A);
       if(solver.info()!=Eigen::Success) {
         // decomposition failed
@@ -400,7 +400,7 @@ namespace Heat_method_3 {
     }
 
 
-    Eigen::VectorXd solve_phi(Matrix c, Matrix divergence, int dimension) const
+    Eigen::VectorXd solve_phi(const Matrix& c, const Matrix& divergence, int dimension) const
     {
 
       Eigen::VectorXd phi;
diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h
index 8347f100084..b3b48587bfd 100644
--- a/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h
+++ b/Heat_method_3/include/CGAL/Heat_method_3/Intrinsic_Delaunay_Triangulation_3.h
@@ -367,6 +367,7 @@ namespace Intrinsic_Delaunay_Triangulation_3 {
            edge_id_map = get(Edge_property_tag(), const_cast<TriangleMesh&>(tm));
            Index edge_i = 0;
            BOOST_FOREACH(edge_descriptor ed, edges(tm)){
+             std::cout << ed << std::endl;
              mark_edges(edge_i,0)=1;
              edge_lengths(edge_i,0) = Polygon_mesh_processing::edge_length(halfedge(ed,tm),tm);
              put(edge_id_map, ed, edge_i++);
@@ -408,7 +409,7 @@ namespace Intrinsic_Delaunay_Triangulation_3 {
              double e2_len = edge_lengths(e2,0);
              double e3_len = edge_lengths(e3,0);
              double angle_a = -(e2_len*e2_len) + e3_len*e3_len + e1_len*e1_len;
-             angle_a = angle_a/(2*e3_len*e1_len);
+             angle_a = acos(angle_a/(2*e3_len*e1_len));
              Point_2 p31(e3_len*std::cos(angle_a), e3_len*std::sin(angle_a));
              put(hcm,hd,p31);
 
diff --git a/Heat_method_3/test/Heat_method_3/heat_method_intrinsic_surface_mesh_test.cpp b/Heat_method_3/test/Heat_method_3/heat_method_intrinsic_surface_mesh_test.cpp
index bbe0c11fa17..fb9244cbcdf 100644
--- a/Heat_method_3/test/Heat_method_3/heat_method_intrinsic_surface_mesh_test.cpp
+++ b/Heat_method_3/test/Heat_method_3/heat_method_intrinsic_surface_mesh_test.cpp
@@ -109,26 +109,28 @@ void check_no_update(const Idt& sm, const Vertex_distance_map& original, const V
 int main(int argc, char*argv[])
 {
   Surface_mesh sm;
-  Vertex_distance_map vdm = get(Vertex_distance_tag(),sm);
-  
-  Idt idt(sm, vdm);
 
-
-  bool idf = false;
-
-  std::ifstream in((argc>1)?argv[1]:"data/pyramid0.off");
+  std::ifstream in((argc>1)?argv[1]:"data/pyramid1.off");
   in >> sm;
   if(!in || num_vertices(sm) == 0) {
     std::cerr << "Problem loading the input data" << std::endl;
     return 1;
   }
 
+  Vertex_distance_map vdm = get(Vertex_distance_tag(),sm);  
+  Idt idt(sm, vdm);
+
+ 
   //source set tests
   Heat_method hm(idt, idt.vertex_distance_map());
   // hm.add_source(* vertices(idt).first);
 
   hm.add_source(boost::graph_traits<Idt>::vertex_descriptor(boost::graph_traits<Surface_mesh>::vertex_descriptor(0),sm));
   hm.update();
+
+  BOOST_FOREACH(boost::graph_traits<Surface_mesh>::vertex_descriptor vd, vertices(sm)){
+    std::cout << get(vdm,vd) << std::endl;
+  }
 #if 0