Update heat_geodesics.cpp (#1512)

* Update heat_geodesics.cpp

* heat geodesic upsampled cube unit test

include/igl/heat_geodesics.cpp

@@ -152,7 +152,7 @@ IGL_INLINE void igl::heat_geodesics_solve(
   }
   const DerivedD div_X = -data.Div*grad_u;
   const DerivedD Beq = (DerivedD(1,1)<<0).finished();
-  igl::min_quad_with_fixed_solve(data.Poisson,(-2.0*div_X).eval(),DerivedD(),Beq,D);
+  igl::min_quad_with_fixed_solve(data.Poisson,(-div_X).eval(),DerivedD(),Beq,D);
   DerivedD Dgamma;
   igl::slice(D,gamma,Dgamma);
   D.array() -= Dgamma.mean();

tests/include/igl/heat_geodesics.cpp

@@ -0,0 +1,27 @@
+#include <test_common.h>
+#include <igl/heat_geodesics.h>
+#include <igl/upsample.h>
+#include <igl/avg_edge_length.h>
+
+TEST_CASE("heat_geodesic: upsampled cube", "[igl]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+  igl::upsample(V,F,3);
+
+  int vid = 0;
+  igl::HeatGeodesicsData<double> data;
+  igl::heat_geodesics_precompute(V,F,data);
+  Eigen::VectorXd dist;
+  igl::heat_geodesics_solve(data, (Eigen::VectorXi(1,1)<<vid).finished(), dist);
+  
+  double avg_edge = igl::avg_edge_length(V,F);
+
+  // Check the adjacent corners
+  for (int i=0; i<6; i++) {
+  REQUIRE((V.row(i)-V.row(0)).norm() == Approx(dist(i)).margin(avg_edge));
+  }
+
+}