diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 6bfbd321b69189614b455189f7df8ccd7ce993ae..7032e02fb52f92760cf26d2b97ea5bc81d0df001 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -688,25 +688,27 @@ void Curvature::computePointareas(){
       const int V0 = _VertexToInt[A->getNum()];
       const int V1 = _VertexToInt[B->getNum()];
       const int V2 = _VertexToInt[C->getNum()];
-      
+
+      factor[0] = 1.0;
+      factor[1] = 1.0;
+      factor[2] = 1.0;
       // if (_isOnBoundary[V0])
       // {
       //     factor[0] = 1.0;
       // }
-      // else { factor[0] = 1.0;}
+      // else {factor[0] = 1.0;}
       // if (_isOnBoundary[V1])
       // {
       //     factor[1] = 1.0;
       // }
       // else {factor[1] = 1.0;}
-
       // if (_isOnBoundary[V2])
       // {
       //     factor[2] = 1.0;
       // }
       // else {factor[2] = 1.0;}
-
-      //Edges
+      
+//Edges
       e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangular element
       e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
       e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
@@ -1128,6 +1130,9 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
       //Estimate curvature based on variations of normals along edges:
       //intialization:
       m = SVector3(0.0, 0.0, 0.0);
+       for (int i = 0; i< 3; ++i)
+	 for (int j = 0; j< 3; ++j)
+	   w(i,j) = 0.0;
 
       //filling:
       for (int j = 0; j< 3; ++j)
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 354eacbdd1ae7d0949dc7a2c43e1f9bfd09f3d02..88dbe9e8c45fab2176f7b751ae51727d17ba32c1 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -1029,7 +1029,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad)
 
 double Centerline::operator() (double x, double y, double z, GEntity *ge)
 {
- 
+
   if (update_needed){
      std::ifstream input;
      input.open(fileName.c_str());
@@ -1041,7 +1041,6 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge)
    }
 
    double xyz[3] = {x,y,z};
-
    //take xyz = closest point on boundary in case we are on the planar in/out faces
    //or in the volume
    bool isCompound = false;
@@ -1078,6 +1077,7 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge)
 
 void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
 {
+
   if (update_needed){
      std::ifstream input;
      input.open(fileName.c_str());
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 7ff2757035e1f5f0db0ad6663b21f41fb0e23089..eae12ec8de278a8bf67ad660e3f731e9ae0002a1 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -369,9 +369,10 @@ void meshGEdge::operator() (GEdge *ge)
     N = ge->meshAttributes.nbPointsTransfinite;
   }
   else{
-    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS || blf)
+    if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || CTX::instance()->mesh.algo2d == ALGO_2D_PACK_PRLGRMS || blf){
       a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
                       CTX::instance()->mesh.lcIntegrationPrecision);
+    }
     else{
       a = Integration(ge, t_begin, t_end, F_Lc, Points,
                       CTX::instance()->mesh.lcIntegrationPrecision);
diff --git a/benchmarks/centerlines/aorta_centerlines.geo b/benchmarks/centerlines/aorta_centerlines.geo
index 694fdceb2d37c523ed41fc0544643dde71e5ffe2..dd0d76ec71b1225a7d2f82393f442e15472de38c 100644
--- a/benchmarks/centerlines/aorta_centerlines.geo
+++ b/benchmarks/centerlines/aorta_centerlines.geo
@@ -1,16 +1,16 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
-Mesh.Algorithm3D = 1; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
+Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
 
-Mesh.LcIntegrationPrecision = 1.e-3;
+Mesh.LcIntegrationPrecision = 1.e-2;
 
-Mesh.RecombineAll = 1;
-Mesh.Bunin = 60;
+//Mesh.RecombineAll = 1;
+//Mesh.Bunin = 60;
 
 Merge "aorta2.stl";
 
 Field[1] = Centerline;
 Field[1].FileName = "centerlinesAORTA.vtk";
-Field[1].nbPoints = 25; //33;
+Field[1].nbPoints = 25; 
 
 Field[1].nbElemLayer = 4;
 Field[1].hLayer = 0.2;//percent of vessel radius
diff --git a/benchmarks/centerlines/carotid_centerlines.geo b/benchmarks/centerlines/carotid_centerlines.geo
index f043b8064495231801b3d4e0eceb4db8128e808b..51aca8bc7a5441a1866ea32b4fd4fb908e3bb73f 100644
--- a/benchmarks/centerlines/carotid_centerlines.geo
+++ b/benchmarks/centerlines/carotid_centerlines.geo
@@ -1,10 +1,10 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 1; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
 Mesh.Algorithm3D = 1;//(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
 
 Mesh.LcIntegrationPrecision = 1.e-2;
 
-Mesh.RecombineAll = 1;
-Mesh.Bunin = 100;
+//Mesh.RecombineAll = 1;
+//Mesh.Bunin = 100;
 
 Merge "carotid.stl";
 
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index 0405dbe15c9b97c02e609805d13d0285d83155ac..8e88efc01c467e184f870e54b97f2ecbd8dac5be 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -1,4 +1,4 @@
-Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
 Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
 
 Mesh.LcIntegrationPrecision = 1.e-2;
@@ -6,7 +6,7 @@ Mesh.LcIntegrationPrecision = 1.e-2;
 Merge "cylemi.stl";
 
 Mesh.RecombineAll = 1;
-//Mesh.Bunin = 200;
+//Mesh.Bunin = 100;
 
 Field[1] = Centerline;
 Field[1].FileName = "centerlinesCYL.vtk";
diff --git a/benchmarks/curvature/test.geo b/benchmarks/curvature/test.geo
index 2e8167cd6dda0f70669c9b4770dd9c825dbb2cc2..d4121badb9a10185409578eadd9ee8abb9bdffbe 100644
--- a/benchmarks/curvature/test.geo
+++ b/benchmarks/curvature/test.geo
@@ -1,4 +1,4 @@
- Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
+ Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
  Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split only with metis 
  
  Mesh.Algorithm=1; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)