diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e8e81abfafb722aeaf3a96140b9ca1956fb36ad2..5c73d488ba43380371d8f20e7625ba6069f0b0e3 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1244,7 +1244,7 @@ void GModel::makeDiscreteRegionsSimplyConnected()
     
     std::vector<std::vector<MElement*> > conRegions;
     int nbRegions = connectedVolumes(allElements, conRegions);
-    remove(*itR);
+    if (nbRegions > 1) remove(*itR);
 
     for(int ire  = 0; ire < nbRegions; ire++){
       int numR = (nbRegions == 1) ? (*itR)->tag() : getMaxElementaryNumber(3) + 1;
@@ -1303,7 +1303,7 @@ void GModel::makeDiscreteFacesSimplyConnected()
     
     std::vector<std::vector<MElement*> > conFaces;
     int nbFaces = connectedSurfaces(allElements, conFaces);
-    remove(*itF);
+    if (nbFaces > 1) remove(*itF);
 
     for(int ifa  = 0; ifa < nbFaces; ifa++){
       int numF = (nbFaces == 1) ? (*itF)->tag() : getMaxElementaryNumber(2) + 1;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 59948c017966a9ddb3cec68dfdfb356cfa9dd852..fb515254a0de658d6424284feb32c59e69592cac 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -475,8 +475,7 @@ GPoint discreteEdge::point(double par) const
 
 SVector3 discreteEdge::firstDer(double par) const
 {
-  
-  double tLoc;
+   double tLoc;
   int iEdge;
   getLocalParameter(par, iEdge, tLoc);
 
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 2d98fdf50ecc396b59e61b92db5034ec78728498..3070425df9475a36cb096cd11c77f6b329898814 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -626,6 +626,31 @@ double gLevelsetQuadric::operator()(const double x, const double y, const double
         + 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
 }
 
+// gLevelsetPopcorn::gLevelsetPopcorn(int tag) : gLevelsetPrimitive(tag) {
+// }
+
+// double gLevelsetPopcorn::operator() (const double x, const double y, const double z) const {
+//   double r0 = 0.25;
+//   double r = sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5));
+//   double  val = r - r0;
+//   for (int k = 0; k< 5; k++){
+//     double xk = r0/sqrt(5)*(2.*cos(2*k*pi/5));
+//     double yk = r0/sqrt(5)*(2.*sin(2*k*pi/5));
+//     double zk = 1.0;
+//     val +=  -4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
+//     }
+//   for (int k = 0; k< 5; k++){
+//     xk = r0/sqrt(5)*(2*cos(2*((k-5)-1)*pi/5));
+//     yk = r0/sqrt(5)*(2*sin(2*((k-5)-1)*pi/5));
+//     zk = 1.0;
+//     val += 4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
+//   }
+//   val  += -4.*exp(((x-0.5-xk)*(x-0.5-xk)+(y-0.5-yk)*(y-0.5-yk)+(z-0.5-zk)*(z-0.5-zk))/0.01);
+//   val  += -4.*math.exp(((x-0.5)**2+(y-0.5)**2+(z-0.5+r0)**2)/0.01);
+//   return val;
+// }
+
+
 gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag) {
     std::vector<std::string> expressions(1, f);
     std::vector<std::string> variables(3);
@@ -657,7 +682,7 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int
 
   double rmax = 0.1;
   double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
-  double rtube = std::max(lx, std::max(ly, lz))*2; //0.5; // * 2.;
+  double rtube = std::max(lx, std::max(ly, lz))*2.;
 
   Msg::Info("Filling coarse point cloud on surfaces");
   std::vector<SPoint3> points;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 63a54f1aec072afdfdc7482aa4dcdb4657fcbe1c..57d58f348d758808f1e81fbb872858148fbaf5f7 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -56,7 +56,7 @@ SMetric3 buildMetricTangentToCurve (SVector3 &t, double curvature, double &lambd
   lambda =  ((2 * M_PI) /( fabs(curvature) *  CTX::instance()->mesh.minCircPoints ) );
 
   SMetric3 curvMetric (1./(lambda*lambda),1.e-10,1.e-10,t,b,c);
-    //SMetric3 curvMetric (1./(lambda*lambda));
+  //SMetric3 curvMetric (1./(lambda*lambda));
   return curvMetric;
 }
 
@@ -233,7 +233,7 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
   case 1:
     {
       GEdge *ged = (GEdge *)ge;
-      Crv = ged->curvature(U)*2;
+      Crv = ged->curvature(U); //*2 WHY A FACTOR 2
       Crv = std::max(Crv, max_surf_curvature(ged, U));
       //Crv = max_surf_curvature(ged, U);      
     }
diff --git a/benchmarks/curvature/Torus/TorusRemesh.geo b/benchmarks/curvature/Torus/TorusRemesh.geo
index 68c370f0e7882061c7872a58293d8639b39478fb..6494f85744f67032ca3f20e10155585c70d44c98 100644
--- a/benchmarks/curvature/Torus/TorusRemesh.geo
+++ b/benchmarks/curvature/Torus/TorusRemesh.geo
@@ -4,15 +4,15 @@
 
 Mesh.Algorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=bamg
 Mesh.RemeshParametrization= 1; //(0) harmonic (1) conformal 
-Mesh.RemeshAlgorithm = 0; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1 
+Mesh.RemeshAlgorithm = 1; //(0) nosplit (1) automatic (2) split only with metis ///Default: 1 
 
 Mesh.CharacteristicLengthFactor=0.4;
 Mesh.CharacteristicLengthFromPoints = 0;
 Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+Mesh.MinimumCirclePoints=45; //default=7
 Mesh.CharacteristicLengthMin = 0.01; //-clmin
-Mesh.CharacteristicLengthMax = 2.0; //-clmax
+Mesh.CharacteristicLengthMax = 0.5; //-clmax
 Mesh.LcIntegrationPrecision = 1.e-5; //-epslc1d
-Mesh.MinimumCirclePoints=15; //default=7
 Mesh.CharacteristicLengthExtendFromBoundary=0;
 
 //Merge "TorusRemeshedBAMG.stl";