diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 22ec69d39422df99b48bc78615127a5f5b20aa23..3b899b5d1a9593d2a0e449f27db065bd8ea7ea49 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2546,6 +2546,7 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
   }
 
   char buffer[256], buffer2[256];
+  std::map<int, std::map<int, std::string> > physicals[4];
 
   if(!fgets(buffer, sizeof(buffer), fp)) return 0; // version line
   if(!fgets(buffer, sizeof(buffer), fp)) return 0; // title
@@ -2695,29 +2696,30 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
     }
   }
   else if (haveLines){
-    printf("have lines \n");
-    int v0, v1;
-    fscanf(fp, "%d", &v0);
-    printf("v0=%d \n", v0);
-    while(fscanf(fp, "%d", &v1) == 1){
-	elements[1][1].push_back(new MLine(vertices[v0],vertices[v1])); 
-    	printf("line =%d %d \n", v0, v1);
-    	v0 = v1;
-    }
-
-    // char str[256] = "XXX";
-    // for (int j = 0; j< numElements; j++){
-    //   printf("line j=%d\n", j);
-    //   if(!fgets(str, sizeof(str), fp)) return 0;
-    //   int v0, v1;
-    //   sscanf(str, "%d", &v0);
-    //   printf("v0=%d \n", v0);
-    //   while(sscanf(str, "%d", &v1) == 1){
-    // 	elements[1][1].push_back(new MLine(vertices[v0],vertices[v1])); 
-    // 	printf("line =%d %d \n", v0, v1);
-    // 	v0 = v1;
-    //   }
-    // }
+    if(!binary){
+      int v0, v1;
+      char line[1024], *p, *pEnd, *pEnd2;
+      int iLine = 1;
+      int num = 0;
+      for (int k= 0; k < numElements; k++){
+	physicals[1][iLine][100] = "centerline";
+	fgets(line, sizeof(line), fp);
+	v0=(int)strtol(line, &pEnd, 10);//ignore firt line
+	v0=(int)strtol(pEnd, &pEnd2, 10);
+	p=pEnd2;
+	while(1){
+	  v1 = strtol(p, &pEnd, 10);
+	  if (p == pEnd )  break;
+	  elements[1][iLine].push_back(new MLine(vertices[v0],vertices[v1])); 
+	  p=pEnd;
+	  v0=v1;
+	}
+	iLine++;
+      }
+    }
+    else{
+      Msg::Error("TODO: implement reading lines for binary files \n");
+    }
   }
 
   for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
@@ -2725,6 +2727,10 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
   _associateEntityWithMeshVertices();
   _storeVerticesInEntities(vertices);
 
+  // store the physical tags
+  for(int i = 0; i < 4; i++)
+    _storePhysicalTagsInEntities(i, physicals[i]);
+
   fclose(fp);
   return 1;
 }
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index cdbde40e8b96e6500133810721fa193f41be13df..dde4c80ff144da3d9416a37c328697afc392e741 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -70,14 +70,6 @@ class MTetrahedron : public MElement {
   {
     return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]);
   }
-  virtual double getMaxEdgeLength(){
-    double maxe = 0.0;
-    for(int i = 0; i < 6; i++){
-      double loce = getEdge(i).length();
-      maxe = std::max(maxe,loce);
-    }
-    return maxe;
-  }
   virtual int getNumEdgesRep(){ return 6; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
   {
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index ddc2525f09f1c8b1e6924661afdc59275988108f..0674c9d137dd8b6047c076e7323d656c25b310b6 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -21,7 +21,6 @@
 #include "BDS.h"
 #include "Context.h"
 #include "GFaceCompound.h"
-#include "SmoothData.h"
 
 #if defined(HAVE_ANN)
 #include "ANN/ANN.h"
@@ -80,7 +79,6 @@ void printVoronoi(GRegion *gr,  std::set<SPoint3> &candidates)
     std::set<MTetrahedron*>::const_iterator it = allTets.begin();
     MTetrahedron *poleTet = *it;
     double maxRadius = poleTet->getCircumRadius();
-    double maxEdgeLength = poleTet->getMaxEdgeLength();
     for(; it != allTets.end(); it++){
       double radius =  (*it)->getCircumRadius();
       if (radius > maxRadius){
@@ -88,7 +86,7 @@ void printVoronoi(GRegion *gr,  std::set<SPoint3> &candidates)
     	poleTet = *it;
       }
     }
-    if (v->onWhat()->dim() == 2 && maxRadius < maxEdgeLength){
+    if (v->onWhat()->dim() == 2 ){
       SPoint3 pc = poleTet->circumcenter();
       // double nx,ny,nz;
       // SVector3 vN = snorm->get(v->x(), v->y(), v->z(), nx,ny,nz);