diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 6e1bd90cbe24ae81d2d23a5cbb591aa17bacb385..cb645f5b64e6fdf54c8ad5145dd3c957b1c9ddb3 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -4,6 +4,7 @@
 # bugs and problems to <gmsh@geuz.org>.
 
 set(SRC
+CenterlineField.cpp
   Generator.cpp
   Field.cpp
     meshGEdge.cpp 
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..88b7791880bf63e72de8d2e56517715b13513ae0
--- /dev/null
+++ b/Mesh/CenterlineField.cpp
@@ -0,0 +1,168 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+//   Emilie Marchandise
+//
+
+#include "CenterlineField.h"
+
+#include <vector>
+#include <map>
+#include <string>
+#include "GModel.h"
+#include "MElement.h"
+#include "MVertex.h"
+#include "MLine.h"
+#include "GEntity.h"
+#include "Field.h"
+
+#if defined(HAVE_ANN)
+#include <ANN/ANN.h>
+
+Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
+
+  index = new ANNidx[1];
+  dist = new ANNdist[1];
+  update_needed = false;
+
+  importFile(fileName);
+  buildKdTree();
+
+
+}
+Centerline::Centerline(): kdtree(0), nodes(0){
+
+  index = new ANNidx[1];
+  dist = new ANNdist[1];
+  update_needed = false;
+
+  options["FileName"] = new FieldOptionString (fileName, "File Name for the centerlines");
+  fileName = "centerlines.vtk";//default
+  //TODO get fileName from file 
+
+  importFile(fileName);
+  buildKdTree();
+
+}
+
+Centerline::~Centerline(){
+  if (mod) delete mod;
+  if(kdtree) delete kdtree;
+  if(nodes) annDeallocPts(nodes);
+  delete[]index;
+  delete[]dist;
+}
+
+void Centerline::importFile(std::string fileName){
+
+  GModel *current = GModel::current();
+
+  mod = new GModel();
+  mod->readVTK(fileName.c_str());
+  mod->removeDuplicateMeshVertices(1.e-8);
+  
+  std::vector<GEntity*> entities ;
+  mod->getEntities(entities) ; 
+    
+  for(unsigned int i = 0; i < entities.size(); i++){
+     if( entities[i]->dim() == 1){
+     for(unsigned int ele = 0; ele < entities[i]->getNumMeshElements(); ele++){ 
+       MLine *l = (MLine*) entities[i]->getMeshElement(ele);
+       MVertex *v0 = l->getVertex(0);
+       MVertex *v1 = l->getVertex(1);
+       std::map<MVertex*, int>::iterator it0 = color.find(v0);
+       std::map<MVertex*, int>::iterator it1 = color.find(v1);
+       if (it0 == color.end() || it1 == color.end()) lines.push_back(l);
+       if (it0 == color.end()) color.insert(std::make_pair(v0, entities[i]->tag()));
+       if (it1 == color.end()) color.insert(std::make_pair(v1, entities[i]->tag()));
+     }
+   }
+ }
+  current->setAsCurrent();
+
+}
+void Centerline::buildKdTree(){
+
+  FILE * f = fopen("myPOINTS.pos","w");
+  fprintf(f, "View \"\"{\n");
+
+  int nbPL = 10;
+  int nbNodes  = (lines.size()+1) + (nbPL*lines.size());
+
+  nodes = annAllocPts(nbNodes, 3);
+  int ind = 0;
+  std::map<MVertex*,int >::iterator itmap = color.begin();
+  while (itmap != color.end()){
+     MVertex *v = itmap->first;
+     nodes[ind][0] = v->x(); 
+     nodes[ind][1] = v->y(); 
+     nodes[ind][2] = v->z(); 
+     itmap++; ind++;
+  }
+  //10 points per line
+ for(unsigned int k = 0; k < lines.size(); ++k){
+   MVertex *v0 = lines[k]->getVertex(0);
+   MVertex *v1 = lines[k]->getVertex(1);
+   SVector3 P0(v0->x(),v0->y(), v0->z());
+   SVector3 P1(v1->x(),v1->y(), v1->z());
+   for (int j= 1; j < nbPL+1; j++){
+     double inc = (double)j/(double)(nbPL+1);
+     SVector3 Pj = P0+inc*(P1-P0);
+     nodes[ind][0] = Pj.x(); 
+     nodes[ind][1] = Pj.y();  
+     nodes[ind][2] = Pj.z();  
+     ind++;
+   }
+ }
+
+ kdtree = new ANNkd_tree(nodes, nbNodes, 3);
+
+ for(unsigned int i = 0; i < nbNodes; ++i){
+   fprintf(f, "SP(%g,%g,%g){%g};\n",
+	   nodes[i][0], nodes[i][1],nodes[i][2],1.0);
+ }
+ fprintf(f,"};\n");
+  fclose(f);
+
+}
+
+double Centerline::operator() (double x, double y, double z, GEntity *ge){
+
+  double xyz[3] = {x,y,z };
+  int num_neighbours = 1;
+  kdtree->annkSearch(xyz, num_neighbours, index, dist);
+  double d = sqrt(dist[0]);
+  double lc = 2*M_PI*d/30.0; //30 mesh elements along circle
+  return lc;
+
+}
+
+void  Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
+  return;
+}
+
+void Centerline::printSplit() const{
+
+  FILE * f = fopen("centerlineSPLIT.pos","w");
+  fprintf(f, "View \"\"{\n");
+
+  for(unsigned int i = 0; i < lines.size(); ++i){
+      MLine *l = lines[i];
+      std::map<MVertex*,int>::const_iterator it0 = 
+        color.find(l->getVertex(0));
+      std::map<MVertex*,int>::const_iterator it1 = 
+        color.find(l->getVertex(1));
+      fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+              l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
+              l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
+  	      (double)it0->second,(double)it1->second);
+   }
+  fprintf(f,"};\n");
+  fclose(f);
+ 
+}
+
+#endif
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
new file mode 100644
index 0000000000000000000000000000000000000000..b8c8ba408991ba96dfaf68fa01d701e87fc4d25a
--- /dev/null
+++ b/Mesh/CenterlineField.h
@@ -0,0 +1,66 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+//   Emilie Marchandise
+
+#ifndef _CENTERLINEFIELD_H_
+#define _CENTERLINEFIELD_H_
+
+#include <vector>
+#include <map>
+#include <string>
+#include "Field.h"
+class GModel;
+class MLine;
+class MVertex;
+class GEntity;
+
+#if defined(HAVE_ANN)
+#include <ANN/ANN.h>
+class ANNkd_tree;
+
+class Centerline : public Field{
+
+  ANNkd_tree *kdtree; 
+  ANNpointArray nodes;
+  ANNidxArray index;
+  ANNdistArray dist;
+  GModel *mod;
+  std::string fileName;
+
+ public:
+  Centerline(std::string fileName);
+  Centerline();
+  ~Centerline();
+
+  void importFile(std::string fileName);
+
+  virtual bool isotropic () const {return false;}
+  virtual const char *getName()
+  {
+    return "centerline Field";
+  }
+  virtual std::string getDescription()
+  {
+    return "The value of this field is the distance to the centerline.\n\n"
+" You should specify a fileName that contains the centerline.\n"
+" The centerline of a surface can be obtained using the open source software vmtk\n\n"
+" using the following script:\n\n"
+"vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
+  }
+  double operator() (double x, double y, double z, GEntity *ge=0);
+  void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
+
+  void printSplit() const;
+  void buildKdTree();
+
+  std::vector<MLine*> lines;
+  std::map<MVertex*,int> color;
+
+};
+#endif
+
+#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 019cbf1a2cbd5e766d651acbfc30ef7534553351..29103eef72985603736ee1f033227759fdd0fddc 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -23,6 +23,7 @@
 #include "Numeric.h"
 #include "mathEvaluator.h"
 #include "BackgroundMesh.h"
+#include "CenterlineField.h"
 
 #if defined(HAVE_POST)
 #include "OctreePost.h"
@@ -92,23 +93,6 @@ class FieldOptionList : public FieldOption
   }
 };
 
-class FieldOptionString : public FieldOption
-{
- public:
-  std::string & val;
-  virtual FieldOptionType getType(){ return FIELD_OPTION_STRING; }
-  FieldOptionString(std::string &_val, std::string _help, bool *_status=0)
-    : FieldOption(_help, _status), val(_val) {}
-  std::string &string() { modified(); return val; }
-  const std::string &string() const { return val; }
-  void getTextRepresentation(std::string &v_str)
-  {
-    std::ostringstream sstream;
-    sstream << "\"" << val << "\"";
-    v_str = sstream.str();
-  }
-};
-
 class FieldOptionPath : public FieldOptionString
 {
  public:
@@ -1931,6 +1915,7 @@ FieldManager::FieldManager()
   map_type_name["Threshold"] = new FieldFactoryT<ThresholdField>();
 #if defined(HAVE_ANN)
   map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
+  map_type_name["Centerline"] = new FieldFactoryT<Centerline>();
 #endif
   map_type_name["Box"] = new FieldFactoryT<BoxField>();
   map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
diff --git a/Mesh/Field.h b/Mesh/Field.h
index c27bd596811f0cb7bd47abf9f3efc343629529be..a21458581e846ec6cbeefd1b6a8cbb1858e33314 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -12,6 +12,11 @@
 #include "GmshConfig.h"
 #include "STensor3.h"
 
+#include <fstream>
+#include <string>
+#include <string.h>
+#include <sstream>
+
 #if defined(HAVE_POST)
 #include "PView.h"
 #endif
@@ -128,4 +133,21 @@ class BoundaryLayerField : public Field {
 };
 #endif
 
+class FieldOptionString : public FieldOption
+{
+ public:
+  std::string & val;
+  virtual FieldOptionType getType(){ return FIELD_OPTION_STRING; }
+  FieldOptionString(std::string &_val, std::string _help, bool *_status=0)
+    : FieldOption(_help, _status), val(_val) {}
+  std::string &string() { modified(); return val; }
+  const std::string &string() const { return val; }
+  void getTextRepresentation(std::string &v_str)
+  {
+    std::ostringstream sstream;
+    sstream << "\"" << val << "\"";
+    v_str = sstream.str();
+  }
+};
+
 #endif
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index 8b5cdf1beb5038aedf176fd31eaca828a1ead229..ec5ef07b648e5a119a89f0ccaf4ea4857b7b88c3 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -105,7 +105,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
     Msg::Error("SLEPc diverged after %d iterations", its);
   else if(reason == EPS_DIVERGED_BREAKDOWN)
     Msg::Error("SLEPc generic breakdown in method");
-#if !(PETSC_VERSION_RELEASE == 0) // petsc-dev
+#if !(PETSC_VERSION_RELEASE == 0 || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 2))) // petsc-dev
   else if(reason == EPS_DIVERGED_NONSYMMETRIC)
     Msg::Error("The operator is nonsymmetric");
 #endif
@@ -154,7 +154,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
       }
        _eigenVectors.push_back(ev);
     }
-#if (PETSC_VERSION_RELEASE == 0) // petsc-dev
+#if (PETSC_VERSION_RELEASE == 0 || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 2))) // petsc-dev
     _try(VecDestroy(&xr));
     _try(VecDestroy(&xi));
 #else
@@ -163,7 +163,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
 #endif
   }
 
-#if (PETSC_VERSION_RELEASE == 0) // petsc-dev
+#if (PETSC_VERSION_RELEASE == 0 || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR == 2))) // petsc-dev
   _try(EPSDestroy(&eps));
 #else
   _try(EPSDestroy(eps));
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 8fe521621f75bf53dddafef6bf011bd9f2c99088..11d3735e1122cefc7eb7523919d198e9004e7554 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -910,8 +910,8 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
  
   //Split the mesh in left and right
   //or Cut the mesh in left and right
-  splitElems(elements); 
-  //cutElems(elements);
+  //splitElems(elements); 
+  cutElems(elements);
 
 }
 
diff --git a/benchmarks/curvature/Torus/TorusInput.stl b/benchmarks/curvature/Torus/TorusInput.stl
index 8fbda00f299c22bace0a0a6b1f7995f2bec4b44e..ec4fa86dc33e5ff9cc6857f27dd89fa4e91d8099 100644
Binary files a/benchmarks/curvature/Torus/TorusInput.stl and b/benchmarks/curvature/Torus/TorusInput.stl differ
diff --git a/gmshpy/gmshMesh.i b/gmshpy/gmshMesh.i
index 7e57df9349ef047230c21fc5912bb43380e48290..cbe71480cfb645b0a2e0fc496d67ab5e82763e52 100644
--- a/gmshpy/gmshMesh.i
+++ b/gmshpy/gmshMesh.i
@@ -13,6 +13,7 @@
   #include "Levy3D.h"
   #include "meshPartition.h"
   #include "meshMetric.h"
+  #include "CenterlineField.h"
 %}
 
 %include std_vector.i
@@ -29,3 +30,4 @@ namespace std {
 %include "Levy3D.h"
 %include "meshPartition.h"
 %include "meshMetric.h"
+%include "CenterlineField.h"