diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index f0d8d20ec5d5c81e093579dc7bd5301e91b52655..11635073a26878a66544e1b31fce2cb2d2053a2b 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -1014,7 +1014,7 @@ void Curvature::writeToPosFile( const std::string & filename)
       //Here is printing the triplet X-Y-Z of each vertex:
       //*************************************************
 
-      outfile << idxelem << " ST("; //VT = vector triangles   //ST = scalar triangle
+      outfile << "ST("; //VT = vector triangles   //ST = scalar triangle
       outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
       outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
       outfile << C->x() << ","<< C->y() << "," << C->z();
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e5f07e70c393b1a5c1620acd44618dbfad05990b..b6f7b69c36ae8acff9a890548be70d06044f49c0 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1394,15 +1394,15 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const
     Curvature& curvature = Curvature::getInstance();
 
     if( !Curvature::valueAlreadyComputed() ) {
-      std::cout << "Need to compute discrete curvature" << std::endl;
-      std::cout << "Getting instance of curvature" << std::endl;
+      Msg::Info("Need to compute discrete curvature for isotropic remesh");
+      Msg::Info("Getting instance of curvature");
 
       curvature.setGModel( model() );
       int computeMax = 0;
       curvature.computeCurvature_Rusinkiewicz(computeMax);
       curvature.writeToPosFile("curvature.pos");
       curvature.writeToVtkFile("curvature.vtk");
-      std::cout << " ... computing curvature finished" << std::endl;
+      Msg::Info(" ... computing curvature finished");
     }
 
     double c0;
@@ -1445,8 +1445,8 @@ double GFaceCompound::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector
     Curvature& curvature = Curvature::getInstance();
 
     if( !Curvature::valueAlreadyComputed() ) {
-      std::cout << "Need to compute discrete curvature" << std::endl;
-      std::cout << "Getting instance of curvature" << std::endl;
+      Msg::Info("Need to compute discrete curvature for anisotropic remesh");
+      Msg::Info("Getting instance of curvature");
 
       curvature.setGModel( model() );
       int computeMax = 0;
@@ -1454,10 +1454,10 @@ double GFaceCompound::curvatures(const SPoint2 &param, SVector3 *dirMax, SVector
       curvature.writeToPosFile("curvature.pos");
       curvature.writeToVtkFile("curvature.vtk");
       curvature.writeDirectionsToPosFile("curvature_directions.pos");
-      std::cout << " ... computing curvature finished" << std::endl;
+      Msg::Info(" ... computing curvature finished");
     }
 
-    std::cout << "I'm using curvatures in GFaceCompound.cpp" << std::endl;
+    //std::cout << "I'm using curvatures in GFaceCompound.cpp" << std::endl;
     double cMin[3];
     double cMax[3];
     SVector3 dMin[3];
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index e484541ac03e23439f87578f2bf1ad6616fc8539..e6d2c612db2bba7fd7b0cd659868473a9867e512 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -509,6 +509,7 @@ double discreteEdge::curvature(double par) const{
     int computeMax = 0;
     curvature.computeCurvature_Rusinkiewicz(computeMax);
     curvature.writeToPosFile("curvature.pos");
+    curvature.writeDirectionsToPosFile("curvature_directions.pos");
     curvature.writeToVtkFile("curvature.vtk");
 
     std::cout << " ... computing curvature finished" << std::endl;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 4db52ed919f87824366af44da0fbf83c60589ae4..1e5b90125df63ff8073f4f926d3d725b744c18bb 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -177,7 +177,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
   std::list<GFace *>::iterator it = faces.begin();
   int count = 0;
   while(it != faces.end()){
-    if ((*it)->geomType() != GEntity::CompoundSurface){
+    if ( ((*it)->geomType() != GEntity::CompoundSurface) && ((*it)->geomType() != GEntity::DiscreteSurface) ){
       SPoint2 par = ge->reparamOnFace((*it), u, 1);
       SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y());
       if (!count) mesh_size = m;
diff --git a/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo b/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
index c0109b161f13ed88fc18cec5af1f456531aa5f77..7943fe5eb7d8b4fd86a357a922c6731b61f04f44 100644
--- a/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
+++ b/benchmarks/curvature/Bifurcation/BifurcationRemeshCurvature.geo
@@ -14,7 +14,7 @@ CreateTopology;
 
 Compound Surface(20) = {1};
 Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal 
-Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1 
+Mesh.RemeshAlgorithm=0; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1 
 
 Physical Surface("Wall") = {20};
 
diff --git a/benchmarks/curvature/Torus/TorusRemesh.geo b/benchmarks/curvature/Torus/TorusRemesh.geo
index 1175e9b522a3de84771b86faddf3e32713c092b7..2e2cb5cfe3794c720248b3711ae4810e48b44ab2 100644
--- a/benchmarks/curvature/Torus/TorusRemesh.geo
+++ b/benchmarks/curvature/Torus/TorusRemesh.geo
@@ -4,7 +4,7 @@
 
 Mesh.Algorithm=7; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=bamg
 Mesh.RemeshParametrization= 1; //(0) harmonic (1) conformal 
-Mesh.RemeshAlgorithm = 1;
+Mesh.RemeshAlgorithm = 0; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1 
 
 Mesh.CharacteristicLengthFactor=0.4;
 Mesh.CharacteristicLengthFromPoints = 0;