diff --git a/Geo/GFace.h b/Geo/GFace.h
index d7d5d559325f83a3b435f35dd4c44a23436732e1..fcaf029cb9b48bef250c8b0a360f6a45bf50c11c 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -338,9 +338,10 @@ class GFace : public GEntity
   // get the boundary layer columns
   BoundaryLayerColumns *getColumns () {return &_columns;}
 
-  std::vector<SPoint3> storage1; //directions storage
-  std::vector<SVector3> storage2; //directions storage
-  std::vector<SVector3> storage3; //directions storage
+  std::vector<SPoint3> storage1; //sizes and directions storage
+  std::vector<SVector3> storage2; //sizes and directions storage
+  std::vector<SVector3> storage3; //sizes and directions storage
+  std::vector<double> storage4; //sizes and directions storage
 };
 
 #endif
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index a782f75c871138f1595bf6a9d1983260bd8d1dcf..97c349256e3f5345d153ec469fe9d2c1333255cc 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -1013,60 +1013,88 @@ void Frame_field::save_energy(GRegion* gr, const std::string& filename){
 Size_field::Size_field(){}
 
 void Size_field::init_region(GRegion* gr){
-  size_t i;
+#if defined(HAVE_ANN)
+  unsigned int i;
   int j;
-  double local_size;
-  double average_x,average_y;
-  SPoint2 point;
+  int index;
+  double h;
+  double e;
+  SPoint3 point;
   MElement* element;
   MVertex* vertex;
   GFace* gf;
   GModel* model = GModel::current();
   std::list<GFace*> faces;
   std::list<GFace*>::iterator it;
-
+  ANNpoint query;
+  ANNidxArray indices;
+  ANNdistArray distances;
+		
   faces = gr->faces();
-
-  boundary.clear();
+	
+  field.clear();
 
   for(it=faces.begin();it!=faces.end();it++){
     gf = *it;
-	if(gf->geomType()==GEntity::CompoundSurface){
-	  ((GFaceCompound*)gf)->deleteInternals();
-	  ((GFaceCompound*)gf)->parametrize();
-	}
-    backgroundMesh::set(gf);
-
-    for(i=0;i<gf->getNumMeshElements();i++){
-      element = gf->getMeshElement(i);
-
-      average_x = 0.0;
-      average_y = 0.0;
 
-      for(j=0;j<element->getNumVertices();j++){
-        vertex = element->getVertex(j);
-        reparamMeshVertexOnFace(vertex,gf,point);
-        average_x = average_x + point.x();
-        average_y = average_y + point.y();
-      }
+	for(i=0;i<gf->storage1.size();i++){
+	  point = gf->storage1[i];
+	  h = gf->storage4[i];
 
-      average_x = average_x/element->getNumVertices();
-      average_y = average_y/element->getNumVertices();
+	  field.push_back(std::pair<SPoint3,double>(point,h));
+	}
+  }
 
-      for(j=0;j<element->getNumVertices();j++){
-        vertex = element->getVertex(j);
-        reparamMeshVertexOnFace(vertex,gf,point);
-        local_size = backgroundMesh::current()->operator()(point.x(),point.y(),0.0);
-        boundary.insert(std::pair<MVertex*,double>(vertex,local_size));
-      }
-    }
+  ANNpointArray duplicate = annAllocPts(field.size(),3);
+	
+  for(i=0;i<field.size();i++){
+    duplicate[i][0] = field[i].first.x();
+	duplicate[i][1] = field[i].first.y();
+	duplicate[i][2] = field[i].first.z();
   }
+	
+  kd_tree = new ANNkd_tree(duplicate,field.size(),3);
 
+  boundary.clear();
+	
+  query = annAllocPt(3);
+  indices = new ANNidx[1];
+  distances = new ANNdist[1];
+  
+  index = 0;
+  e = 0.0;
+	
+  for(it=faces.begin();it!=faces.end();it++){
+    gf = *it;
+    
+	for(i=0;i<gf->getNumMeshElements();i++){
+      element = gf->getMeshElement(i);  
+	
+	  for(j=0;j<element->getNumVertices();j++){
+	    vertex = element->getVertex(j);
+	  
+	    query[0] = vertex->x();
+	    query[1] = vertex->y();
+	    query[2] = vertex->z();
+	
+	    kd_tree->annkSearch(query,1,indices,distances,e);
+	    index = indices[0];
+	
+	    boundary.insert(std::pair<MVertex*,double>(vertex,field[index].second));
+	  }
+	}
+  }
+	
   octree = new MElementOctree(model);
+	
+  annDeallocPt(query);
+  delete[] indices;
+  delete[] distances;
+#endif
 }
 
 void Size_field::solve(GRegion* gr){
-  #if defined(HAVE_PETSC)
+#if defined(HAVE_PETSC)
   linearSystem<double>* system = 0;
   system = new linearSystemPETSc<double>;
 
@@ -1132,7 +1160,7 @@ void Size_field::solve(GRegion* gr){
   }
 
   delete system;
-  #endif
+#endif
 }
 
 double Size_field::search(double x,double y,double z){
@@ -1174,18 +1202,6 @@ double Size_field::search(double x,double y,double z){
   return val;
 }
 
-double Size_field::get_ratio(GFace* gf,SPoint2 point){
-  double val;
-  double uv[2];
-  double tab[3];
-
-  uv[0] = point.x();
-  uv[1] = point.y();
-  buildMetric(gf,uv,tab);
-  val = 1.0/pow(tab[0]*tab[2]-tab[1]*tab[1],0.25);
-  return val;
-}
-
 void Size_field::print_field(GRegion* gr){
   size_t i;
   double min,max;
@@ -1280,9 +1296,15 @@ GRegion* Size_field::test(){
 }
 
 void Size_field::clear(){
-  backgroundMesh::unset();
   delete octree;
+  field.clear();
   boundary.clear();
+#if defined(HAVE_ANN)
+  delete kd_tree->thePoints();
+  delete kd_tree;
+  annClose();
+#endif
+	
 }
 
 /****************class Nearest_point****************/
@@ -1552,7 +1574,7 @@ void Nearest_point::print_field(GRegion* gr){
 	  val = search(x,y,z,vec);
 	  if(val){
 	    print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),
-			  SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
+					  SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
 	  }
 	}
   }
@@ -1579,11 +1601,11 @@ GRegion* Nearest_point::test(){
 void Nearest_point::clear(){
   field.clear();
   vicinity.clear();
-  #if defined(HAVE_ANN)
+#if defined(HAVE_ANN)
   delete kd_tree->thePoints();
   delete kd_tree;
   annClose();
-  #endif
+#endif
 }
 
 /****************static declarations****************/
@@ -1599,8 +1621,12 @@ ANNkd_tree* Frame_field::kd_tree;
 ANNkd_tree* Frame_field::annTree;
 #endif
 
+std::vector<std::pair<SPoint3,double> > Size_field::field;
 std::map<MVertex*,double> Size_field::boundary;
 MElementOctree* Size_field::octree;
+#if defined(HAVE_ANN)
+ANNkd_tree* Size_field::kd_tree;
+#endif
 
 std::vector<SPoint3> Nearest_point::field;
 std::vector<MElement*> Nearest_point::vicinity;
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index 4f7d55b46dd4448c3e3b25adcbefb219fd2cb15a..a89bb17f9ea5d239967ab3b3c16b59ab38ea2221 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -70,14 +70,17 @@ class Frame_field{
 
 class Size_field{
  private:
+  static std::vector<std::pair<SPoint3,double> > field;
   static std::map<MVertex*,double> boundary;
   static MElementOctree* octree;
+#if defined(HAVE_ANN)
+  static ANNkd_tree* kd_tree;
+#endif
   Size_field();
  public:
   static void init_region(GRegion*);
   static void solve(GRegion*);
   static double search(double,double,double);
-  static double get_ratio(GFace*,SPoint2);
   static void print_field(GRegion*);
   static GRegion* test();
   static void clear();
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 5bca57be8e494e62bb91205e380110e2c232cc9d..2fad6b875981f23e23b6576d8b59ea50cf0fbe5c 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2509,6 +2509,7 @@ void directions_storage(GFace* gf){
   int j;
   MVertex* vertex;
   MElement* element;
+  SPoint2 point;
   SVector3 v1;
   SVector3 v2;
   MElementOctree* octree;
@@ -2531,6 +2532,7 @@ void directions_storage(GFace* gf){
   gf->storage1.clear();
   gf->storage2.clear();
   gf->storage3.clear();
+  gf->storage4.clear();
 	
   for(it=vertices.begin();it!=vertices.end();it++){
     ok = 0;
@@ -2548,6 +2550,8 @@ void directions_storage(GFace* gf){
       gf->storage1.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
 	  gf->storage2.push_back(v1);
 	  gf->storage3.push_back(v2);
+	  reparamMeshVertexOnFace(*it,gf,point);
+	  gf->storage4.push_back(backgroundMesh::current()->operator()(point.x(),point.y(),0.0));
 	}  
   }