diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index b6f7b69c36ae8acff9a890548be70d06044f49c0..8922a80edaa5b338c6861cc6ce27dec1055c9566 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -83,6 +83,19 @@ int intersection_segments (SPoint3 &p1, SPoint3 &p2,
 	    x[1] >= 0.0 && x[1] <= 1.);
   }
 }
+static void printBound(std::vector<MVertex*> &l, int tag){
+  //print boundary
+  char name[256];
+  sprintf(name, "myBOUND%d.pos", tag);
+  FILE * xyz = fopen(name,"w");
+  fprintf(xyz,"View \"\"{\n");
+  for(int i = 0; i < l.size(); i++){
+   MVertex *v = l[i];
+   fprintf(xyz,"SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), i);
+ }
+ fprintf(xyz,"};\n");
+ fclose(xyz);
+}
 
 static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
                           std::vector<double> &coord)
@@ -144,6 +157,7 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
     }
     if(!found) return false;
   }    
+
   return true;
 }
 
@@ -628,7 +642,8 @@ bool GFaceCompound::parametrize() const
   if(allNodes.empty()) buildAllNodes();
   
   bool success = orderVertices(_U0, _ordered, _coords);
-  if(!success) Msg::Error("Could not order vertices on boundary");
+  printBound(_ordered, this->tag());
+  if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
     
   // Laplace parametrization
   if (_mapping == HARMONIC){
@@ -636,6 +651,8 @@ bool GFaceCompound::parametrize() const
     fillNeumannBCS();
     parametrize(ITERU,HARMONIC); 
     parametrize(ITERV,HARMONIC);
+    //parametrize(ITERU,CONVEXCOMBINATION);
+    //parametrize(ITERV,CONVEXCOMBINATION);
   }
   // Multiscale Laplace parametrization
   else if (_mapping == MULTISCALE){
@@ -663,9 +680,9 @@ bool GFaceCompound::parametrize() const
   // Radial-Basis Function parametrization
   else if (_mapping == RBF){    
     Msg::Debug("Parametrizing surface %d with 'RBF' ", tag());
-     
+    double t1 = Cpu();
     int variableEps = 0;
-    int radFunInd = 1; //MQ RBF
+    int radFunInd = 1; // 1 MQ RBF , 0 GA
     double sizeBox = getSizeH();
 
     fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
@@ -683,7 +700,8 @@ bool GFaceCompound::parametrize() const
     //_rbf->computeCurvature(coordinates);
     //printStuff();
     //exit(1);
-
+    double t2 = Cpu();
+    printf("param RBF in %g sec \n", t2-t1);
   }
 
   buildOct();  
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 94195882f6738208ea46f4d00dc2d9643693695d..52f9610d754874114d110be37fa78adfc5c4bffc 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -74,13 +74,13 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
 {
 
   allCenters.resize(allNodes.size(),3);
-  double tol =  4.e-2*sBox;
+  double tol =  6.e-2*sBox;
 
   //remove duplicate vertices
   //add bc nodes
   for(unsigned int i = 0; i < bcNodes.size(); i++){
     myNodes.insert(bcNodes[i]);
-    //if (bcNodes.size()  > 20) i+=3;
+    //if (bcNodes.size()  > 20) i+=2;
   }
   
   //then  create Mvertex position
@@ -104,6 +104,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
   normals.resize(nbNodes,3);
   int index = 0;
   double dist_min = 1.e6;
+  double dist_max = 1.e-6;
   for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end() ; ++itv){
     MVertex *v1 = *itv;
     centers(index,0) = v1->x()/sBox;
@@ -118,12 +119,12 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
       MVertex *v2 = *itv2;
       double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+(v1->y()-v2->y())*(v1->y()-v2->y())
 			 +(v1->z()-v2->z())*(v1->z()-v2->z()))/sBox;
-      if (dist<dist_min && *itv != *itv2 && dist > 1.e-6) dist_min = dist;
+      if (dist<dist_min && *itv != *itv2 && dist > 1.e-5) dist_min = dist;
     }
     index++;
   }
  
-  delta = 0.33*dist_min;//0.6
+  delta = 0.33*dist_min;//0.33
   radius= nbNodes/10.;   
 
   matAInv.resize(nbNodes, nbNodes);
@@ -134,7 +135,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
 
   Msg::Info("*****************************************");
   Msg::Info("*** RBF nbNodes=%d (all=%d) ", myNodes.size(), allNodes.size());
-  Msg::Info("*** RBF rad=%g dist_min =%g ", radius, dist_min);
+  Msg::Info("*** RBF rad=%g dist_min =%g", radius, dist_min);
   Msg::Info("*****************************************");  
 
   printNodes(myNodes);
@@ -843,7 +844,7 @@ bool GRbf::UVStoXYZ(const double  u_eval, const double v_eval,
     for (int j = i+1; j < num_neighbours; j++){
       double dist = sqrt((UV(index[i],0)-UV(index[j],0))*(UV(index[i],0)-UV(index[j],0))+
     			 (UV(index[i],1)-UV(index[j],1))*(UV(index[i],1)-UV(index[j],1)));
-      if (dist<dist_min && dist > 1.e-6) dist_min = dist;
+      if (dist<dist_min && dist > 1.e-5) dist_min = dist;
     }
  }
  delete [] index;
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 5f9f3d9a71d2d78f6f019b5a031244b93869c554..f589a260a73ee71c8f266738a84805e818421307 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -467,7 +467,6 @@ class functionLevelsetSmooth : public function {
   double _valMin, _valPlus, _E;
   void call(dataCacheMap *m, fullMatrix<double> &val) 
   {
-    
     double ivalPlus = 1./_valPlus;
     double ivalMin = 1./_valMin;
     for (int i = 0; i < val.size1(); i++)
@@ -491,6 +490,7 @@ class functionLevelsetSmooth : public function {
   }
   functionLevelsetSmooth(const function *f0, const double valMin, const double valPlus, const double E) : function(f0->getNbCol()) 
   {
+    printf("Levelset bandwidth is E = %g \n", E);
     setArgument (_f0, f0);
     _valMin  = valMin;
     _valPlus = valPlus;