diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1568dcba48937b17a6d7b7fe2d2ea998b96189d8..e88a83ff0c8bf0bcad1c8cc6646edb3d71db20b8 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1183,11 +1183,11 @@ bool GFaceCompound::parametrize_conformal_spectral() const
     myAssembler.numberVertex(v, 0, 2);
   }
   
-  for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
-    MVertex *v = *itv;
-    myAssembler.numberVertex(v, 0, 1);
-    myAssembler.numberVertex(v, 0, 2);
-  }
+  // for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+  //   MVertex *v = *itv;
+  //   myAssembler.numberVertex(v, 0, 1);
+  //   myAssembler.numberVertex(v, 0, 2);
+  // }
   
   laplaceTerm laplace1(model(), 1, ONE);
   laplaceTerm laplace2(model(), 2, ONE);
@@ -1220,13 +1220,13 @@ bool GFaceCompound::parametrize_conformal_spectral() const
       myAssembler.assemble(v, 0, 2, v, 0, 2,  epsilon);
     }
   }
-  for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
-    MVertex *v = *itv;
-    if (std::find(_ordered.begin(), _ordered.end(), v) == _ordered.end() ){
-      myAssembler.assemble(v, 0, 1, v, 0, 1,  epsilon);
-      myAssembler.assemble(v, 0, 2, v, 0, 2,  epsilon);
-    }
-  }
+  // for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+  //   MVertex *v = *itv;
+  //   if (std::find(_ordered.begin(), _ordered.end(), v) == _ordered.end() ){
+  //     myAssembler.assemble(v, 0, 1, v, 0, 1,  epsilon);
+  //     myAssembler.assemble(v, 0, 2, v, 0, 2,  epsilon);
+  //   }
+  // }
   
   myAssembler.setCurrentMatrix("B");
   
@@ -1249,7 +1249,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
   // FIXME: force non-hermitian. For some reason (roundoff errors?)
   // on some machines and for some meshes slepc complains about bad IP
   // norm otherwise
-  eigenSolver eig(&myAssembler, "B" , "A", true);
+  eigenSolver eig(&myAssembler, "B" , "A", false);
   bool converged = eig.solve(1, "largest");
     
   if(converged) {
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index a0b1366f7de79cddc2f7d7c0b495b41f5b3c90cd..edc5a94871bbedcaec75e531d79fa56530dd4322 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -1035,16 +1035,60 @@ void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
    //double lc = operator()(x,y,z,ge);
    //metr = SMetric3(1./(lc*lc));
 
-   // double xyz[3] = {x,y,z };
-   // ANNidxArray index2 = = new ANNidx[1];
-   // ANNdistArray dist2 =  new ANNdist[1];
-   // int num_neighbours = 2;
-   // kdtree->annkSearch(xyz, num_neighbours, index2, dist2);
-   // double d = sqrt(dist2[0]);
-   // MVertex *v = vertices[ind2s[i]];  
-
-   // delete[]index2;
-   // delete[]dist2; 
+   double xyz[3] = {x,y,z };
+   ANNidxArray index2 = new ANNidx[2];
+   ANNdistArray dist2 = new ANNdist[2];
+   int num_neighbours = 2;
+   kdtree->annkSearch(xyz, num_neighbours, index2, dist2);
+   double d = sqrt(dist2[0]);
+   double lc = 2*M_PI*d/nbPoints; 
+   SVector3  p0(nodes[index2[0]][0], nodes[index2[0]][1], nodes[index2[0]][2]);
+   SVector3  p1(nodes[index2[1]][0], nodes[index2[1]][1], nodes[index2[1]][2]);
+   SVector3 dir = p1-p0;
+   dir.normalize();
+   SVector3 dir1, dir2;
+   if (dir[1]!=0.0 && dir[2]!=0.0){
+     dir1 = SVector3(1.0, 0.0, -dir[0]/dir[2]);
+     dir2 = SVector3 (dir[0]/dir[2], -(dir[0]*dir[0]+dir[2]*dir[2])/(dir[1]*dir[2]), 1.0);
+   }
+   else if (dir[0]!=0.0 && dir[2]!=0.0){
+     dir1 = SVector3(-dir[1]/dir[0], 1.0, 0.0);
+     dir2 = SVector3(1.0, dir[1]/dir[0], -(dir[1]*dir[1]+dir[0]*dir[0])/(dir[0]*dir[2]));
+   }
+   else if (dir[0]!=0.0 && dir[1]!=0.0){
+     dir1 = SVector3(0.0, -dir[2]/dir[1], 1.0);
+     dir2 = SVector3(-(dir[1]*dir[1]+dir[2]*dir[2])/(dir[0]*dir[1]), 1.0, dir[2]/dir[1]);
+   }
+   else if (dir[0]==0.0 && dir[1]==0.0){
+     dir1 = SVector3(0.0, 1.0, 0.0);
+     dir2 = SVector3(1.0, 0.0, 0.0);
+   }
+   else if (dir[1]==0.0 && dir[2]==0.0){
+     dir1 = SVector3(0.0, 1.0, 0.0);
+     dir2 = SVector3(0.0, 0.0, 1.0);
+   }
+   else if (dir[0]==0.0 && dir[2]==0.0){
+     dir1 = SVector3(1.0, 0.0, 0.0);
+     dir2 = SVector3(0.0, 0.0, 1.0);
+   }
+   else {printf("ARGHH EMI DO STHG \n"); exit(1);}
+   // printf("XYZ =%g %g %g r=%g dir0 = %g %g %g dir1 = %g %g %g dir2 =%g %g %g\n", 
+   // 	  x,y,z,d, dir[0], dir[1], dir[2], dir1[0], dir1[1], dir1[2],  dir2[0], dir2[1], dir2[2] );
+   // printf("0x1 =%g 1x2=%g 2x1=%g \n", dot(dir, dir1), dot(dir1,dir2), dot(dir2,dir));
+   dir1.normalize();
+   dir2.normalize();
+   
+   //dir = SVector3(1.0, 0.0, 0.0);
+   //dir1 = SVector3(0.0, 1.0, 0.0);
+   //dir2 = SVector3(0.0, 0.0, 1.0);
+   
+   double lcA = 4.*lc;
+   double lam1 = 1./(lcA*lcA);
+   double lam2 = 1./(lc*lc);
+   metr = SMetric3(lam1,lam2,lam2, dir, dir1, dir2);
+
+   delete[]index2;
+   delete[]dist2; 
 
    return;
 }
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 0c06e14be5e1d1ca97c56bac05537857d280a403..876b24657ca2a75bb177e378df0e2d0d5cc74182 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -30,6 +30,7 @@
 #include "meshGFaceLloyd.h"
 #include "CenterlineField.h"
 #include "Field.h"
+#include "Options.h"
 
 #if defined(HAVE_POST)
 #include "PView.h"