diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 42625d1f57298d31f8c15a69de83437cc997aa34..dfa238787b2156cc82a8175e0f066d2aab593db5 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -34,6 +34,9 @@
 #include "OpenFile.h"
 #include "CreateFile.h"
 #include "Options.h"
+#include "meshGEdge.h"
+#include "meshGFace.h"
+#include "meshGRegion.h"
 
 #if defined(HAVE_MESH)
 #include "Field.h"
@@ -516,81 +519,105 @@ int GModel::mesh(int dimension)
 }
 
 
-int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<double> parameters)
+int GModel::adaptMesh(int technique, simpleFunction<double> *f, std::vector<double> parameters, bool meshAll)
 {
 #if defined(HAVE_MESH)
 
-  if (getNumMeshElements() == 0) mesh(getDim());
-  meshMetric *mm; 
- 
   int ITER = 0;
   int nbElemsOld = 0;
   int nbElems;
-  while(1){
-    Msg::Info("-- adaptMesh ITER =%d ", ITER);
-    std::vector<MElement*> elements;
-
-    if (getDim() == 2){
-      for (fiter fit = firstFace(); fit != lastFace(); ++fit){
-	if ((*fit)->quadrangles.size())return -1;
-	for (int i=0;i<(*fit)->triangles.size();i++){
-	  elements.push_back((*fit)->triangles[i]);
+  int niter = parameters.size() >=4 ? (int) parameters[3] : 3;
+
+  if (meshAll){
+  
+    meshMetric *bgm = 0;
+    FieldManager *fields = getFields();
+    fields->reset();
+    while(1){
+      Msg::Info("-- adaptMesh (allDim) ITER =%d ", ITER);
+
+      opt_mesh_algo2d(0, GMSH_SET, 7.0); //bamg
+      opt_mesh_algo3d(0, GMSH_SET, 7.0); //mmg3d
+      GenerateMesh(this, getDim());
+      nbElems = getNumMeshElements();
+
+      if (fields) fields->deleteField(1);
+      if (++ITER >= niter)  break;
+      if (ITER > 5 && fabs((double)(nbElems - nbElemsOld)) < 0.005 * nbElemsOld) break;
+	
+      //if(bgm) delete bgm ; //do not do this since we have already deleted the field
+      bgm = new meshMetric(this, technique, f, parameters);
+      int id = fields->newId();
+      (*fields)[id] = bgm;
+      fields->background_field = id;
+      fields->printField();
+      
+      std::for_each(firstEdge(), lastEdge(), deMeshGEdge());
+      std::for_each(firstFace(), lastFace(), deMeshGFace());
+      std::for_each(firstRegion(), lastRegion(), deMeshGRegion());
+
+      nbElemsOld = nbElems;
+      
+    }
+    //if (bgm) delete bgm;
+  }
+
+  else{
+
+    if (getNumMeshElements() == 0) mesh(getDim());
+    meshMetric *mm; 
+
+    while(1) {
+      Msg::Info("-- adaptMesh ITER =%d ", ITER);
+      std::vector<MElement*> elements;
+
+      if (getDim() == 2){
+	for (fiter fit = firstFace(); fit != lastFace(); ++fit){
+	  if ((*fit)->quadrangles.size())return -1;
+	  for (int i=0;i<(*fit)->triangles.size();i++){
+	    elements.push_back((*fit)->triangles[i]);
+	  }
 	}
       }
-    }
-    else if (getDim() == 3){
-      for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
-	if ((*rit)->hexahedra.size())return -1;
-	for (int i=0;i<(*rit)->tetrahedra.size();i++){
-	  elements.push_back((*rit)->tetrahedra[i]);
+      else if (getDim() == 3){
+	for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
+	  if ((*rit)->hexahedra.size())return -1;
+	  for (int i=0;i<(*rit)->tetrahedra.size();i++){
+	    elements.push_back((*rit)->tetrahedra[i]);
+	  }
 	}
       }
-    }
 
-    nbElems = elements.size();
-    if (nbElems == 0)return -1;
+      nbElems = elements.size();
+      if (nbElems == 0)return -1;
 
-    double lcmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
-    double lcmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
-    int niter = parameters.size() >=4 ? (int) parameters[3] : 3;
+      mm = new meshMetric(this, technique, f, parameters);
+      mm->setAsBackgroundMesh (this);
 
-    switch(technique){
-    case 1 : 
-      mm = new meshMetric (elements, f, meshMetric::LEVELSET,lcmin,lcmax,parameters[0], 0);
-      break;
-    case 2 :
-      mm = new meshMetric (elements, f, meshMetric::HESSIAN,lcmin,lcmax,0, parameters[0]);
-      break;
-    case 3 :
-      mm = new meshMetric (elements, f, meshMetric::FREY,lcmin,lcmax,parameters[0], 0.01);
-      break;
-    default : Msg::Error("Unknown Adaptive Strategy");return -1;
-    }
-    // the background mesh is the mesh metric
-    mm->setAsBackgroundMesh (this);
-    if (getDim() == 2){
-      for (fiter fit = firstFace(); fit != lastFace(); ++fit){
-	if((*fit)->geomType() != GEntity::DiscreteSurface){
-	  meshGFaceBamg(*fit);
-	  //laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
+      if (getDim() == 2){
+	for (fiter fit = firstFace(); fit != lastFace(); ++fit){
+	  if((*fit)->geomType() != GEntity::DiscreteSurface){
+	    meshGFaceBamg(*fit);
+	    laplaceSmoothing(*fit,CTX::instance()->mesh.nbSmoothing);
+	  }
+	  if(_octree) delete _octree;
+	  _octree = 0;
 	}
-	if(_octree) delete _octree;
-	_octree = 0;
       }
-    }
-    else if (getDim() == 3){
-      for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
-	refineMeshMMG(*rit);
-	if(_octree) delete _octree;
-	_octree = 0;
+      else if (getDim() == 3){
+	for (riter rit = firstRegion(); rit != lastRegion(); ++rit){
+	  refineMeshMMG(*rit);
+	  if(_octree) delete _octree;
+	  _octree = 0;
+	}
       }
-    }
-    delete mm;
-    if (++ITER >= niter) break;
-    if (fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
+      delete mm;
+      if (++ITER >= niter) break;
+      if (fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
 
-    nbElemsOld = nbElems;
+      nbElemsOld = nbElems;
 
+    }
   }
 
   return 0;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index adf35614b3d68833d818e306ba4648a462414e61..6bf05e1046fcfe8002c171df132157df964c38b3 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -382,7 +382,7 @@ class GModel
   // we assume that boundaries are already adapted.
   // This should be fixed.
   
-  int adaptMesh (int technique, simpleFunction<double> *f, std::vector<double> parameters);
+  int adaptMesh (int technique, simpleFunction<double> *f, std::vector<double> parameters, bool meshAll=false);
 
   // make the mesh a high order mesh at order N
   // linear is 1 if the high order points are not placed on the geometry of the model 
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index f63625caddb09087ec6740ebb5eeed7bc89faf55..067f99f3ac9095b8b69748eb2a1d8424459b6c5b 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1095,7 +1095,7 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
       if(vertexMap.count(numV))
         vmv.push_back(vertexMap[numV]);
       else {
-        MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV);
+        MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV); 
         vmv.push_back(mv);
         vertexMap[numV] = mv;
       }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index d19e1e6a5f776590957780f55339cb1c36db7d05..e46a6547ebc0b1ee9bf5f5c75d9c01ce2b3b610d 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -340,6 +340,7 @@ class MElement
   virtual MElement *copy(std::map<int, MVertex*> &vertexMap,
                          std::map<MElement*, MElement*> &newParents,
                          std::map<MElement*, MElement*> &newDomains);
+
 };
 
 class MElementFactory{
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 9c96953e866c55c158ff4692fec98780f4e28171..1240ff9e4c8cf5f5395445f134793305037ca604 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -100,6 +100,7 @@ class FieldManager : public std::map<int, Field*> {
   int boundaryLayer_field;
   // compatibility with -bgm
   void setBackgroundMesh(int iView);
+  inline void printField(){ for(std::map<int, Field *>::iterator it = begin(); it != end(); it++) printf("id=%d %d\n", it->first);};
 };
 
 // Boundary Layer Field (used both for anisotropic meshing and BL
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 3d570f78544f5f86600813ffc40eed2784d35502..1bddd8249a1fb89d8896283e75ea907a4d1be233 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1123,7 +1123,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
       bowyerWatson(gf);
     else {
       bowyerWatson(gf);
-      printf("in bamg *** \n");
       meshGFaceBamg(gf);
     }
     laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index eabd8c80fd83d6f5b3017d4b596aacb0c831e56c..67212185ec9eb74784ce92f6f27bac7a64661132 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -76,7 +76,7 @@ void meshGFaceBamg(GFace *gf){
   Vertex2 *bamgVertices = new Vertex2[all.size()];
   int index = 0;
   for(std::set<MVertex*>::iterator it = all.begin(); it!=all.end(); ++it){
-    if ((*it)->onWhat()->dim() < 2){
+    if ((*it)->onWhat()->dim() <= 1){
       SPoint2 p;
       bool success = reparamMeshVertexOnFace(*it, gf, p);
       bamgVertices[index][0] = p.x();
@@ -123,8 +123,8 @@ void meshGFaceBamg(GFace *gf){
   for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
       numEdges += (*it)->lines.size();
   }
-  Seg *bamgBoundary = new Seg[numEdges];
 
+  Seg *bamgBoundary = new Seg[numEdges];
   int count = 0;
   for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
     for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index 4d9533e4de520e347b19fe120fe958fa182ae83e..634676a187fe6317c7dfd734964c2e4a1122a65b 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -271,6 +271,8 @@ void refineMeshMMG(GRegion *gr)
     Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------");
     mmg3d::MMG_mmg3dlib(opt,mmg,sol); 
     Msg::Debug("-------- MG3D TERMINATED -------------------");
+    Msg::Info("MMG3D succeeded %d vertices %d tetrahedra",
+	      mmg->np, mmg->ne);
     // Here we should interact with BGM
     updateSizes(gr,mmg, sol);
 
@@ -278,8 +280,8 @@ void refineMeshMMG(GRegion *gr)
     if (fabs((double)(nTnow - nT)) < 0.05 * nT) break;
   }  
 
-  char test[] = "test.mesh";  
-  MMG_saveMesh(mmg, test);
+  //char test[] = "test.mesh";  
+  //MMG_saveMesh(mmg, test);
 
   gr->deleteVertexArrays();
   for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index ae8ca9dde25ad2cede7fc82d24f3fbbe487b610a..16f11ec54b9102efef9183c25baec09b95798403 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -4,6 +4,7 @@
 #include "GEntity.h"
 #include "GModel.h"
 #include "gmshLevelset.h"
+#include "MElementOctree.h"
 
 static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &lt){
   std::set<MElement*> stencil;
@@ -22,14 +23,47 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
   lt.insert(lt.begin(),stencil.begin(),stencil.end());
 }
 
-meshMetric::meshMetric ( std::vector<MElement*> &elements, simpleFunction<double> *fct,  meshMetric::MetricComputationTechnique technique, double lcmin, double lcmax, double E, double eps)
-  : _technique(technique), _epsilon(eps), _E(E), _fct(fct), hmax(lcmax),hmin(lcmin) 
-{
-  _dim = elements[0]->getDim();
-  computeMetric(elements);
+meshMetric::meshMetric(GModel *gm, int technique, simpleFunction<double> *fct, std::vector<double> parameters): _fct(fct){
+
+  _dim  = gm->getDim();
+  std::map<int, MVertex*> vertexMap;
+  std::map<MElement*, MElement*> newP;
+  std::map<MElement*, MElement*> newD;
+
+  if (_dim == 2){
+    for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); ++fit){
+      for (int i=0;i<(*fit)->getNumMeshElements();i++){
+  	MElement *e = (*fit)->getMeshElement(i);
+  	MElement *copy = e->copy(vertexMap, newP, newD);
+  	_elements.push_back(copy);
+      }
+    }
+  }
+  else if (_dim == 3){
+    for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
+      for (int i=0;i<(*rit)->getNumMeshElements();i++){
+  	MElement *e = (*rit)->getMeshElement(i);
+  	MElement *copy = e->copy(vertexMap, newP, newD);
+  	_elements.push_back(copy);
+      }
+    }
+  }
+
+  _octree = new MElementOctree(_elements); 
+  hmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
+  hmax = parameters.size() >=3 ? parameters[2] : CTX::instance()->mesh.lcMax;
+  
+  _E = parameters[0];
+  _epsilon = parameters[0];
+  _technique =  (MetricComputationTechnique)technique;
+ 
+  computeMetric();
   //printMetric("toto.pos");
-}
 
+}
+meshMetric::~meshMetric(){
+  if (_octree) delete _octree;
+}
 
 void meshMetric::computeValues( v2t_cont adj){
 
@@ -51,12 +85,12 @@ void meshMetric::computeHessian( v2t_cont adj){
     while (it != adj.end()) {
       std::vector<MElement*> lt = it->second;      
       MVertex *ver = it->first;
-      while (lt.size() < 7) increaseStencil(ver,adj,lt);
-      if ( ver->onWhat()->dim() < _dim ){
-	while (lt.size() < 12){
-	  increaseStencil(ver,adj,lt);
-	}
-      }
+      while (lt.size() < 10) increaseStencil(ver,adj,lt);
+      // if ( ver->onWhat()->dim() < _dim ){
+      // 	while (lt.size() < 12){
+      // 	  increaseStencil(ver,adj,lt);
+      // 	}
+      // }
       
       fullMatrix<double> A  (lt.size(),DIM);
       fullMatrix<double> AT (DIM,lt.size());
@@ -103,12 +137,12 @@ void meshMetric::computeHessian( v2t_cont adj){
   }
    
 }
-void meshMetric::computeMetric(std::vector<MElement*> &e){
+void meshMetric::computeMetric(){
   
   v2t_cont adj;
-  buildVertexToElement (e,adj);
+  buildVertexToElement (_elements,adj);
 
-  //printf("%d elements are considered in the metric \n",e.size());
+  printf("%d elements are considered in the meshMetric \n",_elements.size());
 
   computeValues(adj);
   computeHessian(adj);
@@ -142,7 +176,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
        double divEps = 1./0.01;
        double norm = gr(0)*gr(0)+gr(1)*gr(1)+gr(2)*gr(2);
        if (dist < _E && norm != 0.0){
-	 double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
+	 double h = hmin*(hmax/hmin-1.0)*dist/_E + hmin;
 	 double C = 1./(h*h) -1./(hmax*hmax);
 	 /*	 hfrey(0,0) += C*gr(0)*gr(0)/(norm) + hessian(0,0)*divEps; //metric intersection ???
 		 hfrey(1,1) += C*gr(1)*gr(1)/(norm) + hessian(1,1)*divEps;
@@ -154,16 +188,15 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
 	 hfrey(0,0) += C*gr(0)*gr(0)/(norm) ;
 	 hfrey(1,1) += C*gr(1)*gr(1)/(norm) ;
 	 hfrey(2,2) += C*gr(2)*gr(2)/(norm) ;
-	 hfrey(1,0) = C*gr(1)*gr(0)/(norm) ;
-	 hfrey(2,0) = C*gr(2)*gr(0)/(norm) ;
-	 hfrey(2,1) = C*gr(2)*gr(1)/(norm) ;	 	 
+	 hfrey(1,0) = hfrey(0,1) = C*gr(1)*gr(0)/(norm) ;
+	 hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) ;
+	 hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) ;	 	 
        }
        SMetric3 sss; sss.setMat(hessian);
        sss *= divEps;
        sss(0,0) += 1/(hmax*hmax);
        sss(1,1) += 1/(hmax*hmax);
        sss(2,2) += 1/(hmax*hmax);
-
        H = intersection(sss,hfrey);
      }
      //See paper Hachem and Coupez: 
@@ -171,7 +204,7 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
      //industrial furnaces using the immersed volume technique, ijnmf, 2010
      else if (_technique == meshMetric::LEVELSET ){
       SVector3 gr = grads[ver];
-      fullMatrix<double> hlevelset(3,3);
+      fullMatrix<double> hlevelset(3,3); 
       hlevelset(0,0) = 1./(hmax*hmax);
       hlevelset(1,1) = 1./(hmax*hmax);
       hlevelset(2,2) = 1./(hmax*hmax);
@@ -273,12 +306,10 @@ void meshMetric::computeMetric(std::vector<MElement*> &e){
 }
 
 double meshMetric::operator() (double x, double y, double z, GEntity *ge){
-  GModel *gm = _nodalMetrics.begin()->first->onWhat()->model();
-  SPoint3 xyz(x,y,z),uvw;
-
+  SPoint3 xyz(x,y,z), uvw;
   double initialTol = MElement::getTolerance();
   MElement::setTolerance(1.e-4); 
-  MElement *e = gm->getMeshElementByCoord(xyz,_dim);    
+  MElement *e = _octree->find(x, y, z, _dim);  
   MElement::setTolerance(initialTol);
   if (e){
     e->xyz2uvw(xyz,uvw);
@@ -295,14 +326,12 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge){
 
 void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
   metr = SMetric3(1.e-22);
-  GModel *gm = _nodalMetrics.begin()->first->onWhat()->model();
-  SPoint3 xyz(x,y,z),uvw;
+  SPoint3 xyz(x,y,z), uvw;
   double initialTol = MElement::getTolerance();
   MElement::setTolerance(1.e-4); 
-  MElement *e = gm->getMeshElementByCoord(xyz,_dim);    
+  MElement *e = _octree->find(x, y, z, _dim);   
   MElement::setTolerance(initialTol);
 
-
   if (e){
     e->xyz2uvw(xyz,uvw);
     SMetric3 m1 = _nodalMetrics[e->getVertex(0)];
@@ -320,7 +349,6 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
   }
 }
 
-// FIXME !!!!!!!!! 2 be fixed !!!
 void meshMetric::setAsBackgroundMesh (GModel *gm){
   FieldManager *fields = gm->getFields();
   int id = fields->newId();
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 240ad8d9833d599cdbfdc234692dc85ba3a2e323..84b63e5d3c194f0f9280a4a529ea54ab69b1bb99 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -7,11 +7,12 @@
 template <class scalar> class simpleFunction;
 class MVertex;
 class gLevelset;
+class MElementOctree;
 
 /**Anisotropic mesh size field based on a metric */
 class meshMetric: public Field {
  public:
-  typedef enum {FREY=1,LEVELSET=2,HESSIAN=3} MetricComputationTechnique;
+  typedef enum {LEVELSET=1,HESSIAN=2, FREY=3} MetricComputationTechnique;
  private:
   int _dim;
   double _epsilon, _E;
@@ -19,6 +20,9 @@ class meshMetric: public Field {
   double hmin, hmax;
   simpleFunction<double> *_fct;
 
+  std::vector<MElement*> _elements;
+  MElementOctree *_octree;
+
   std::map<MVertex*,double> vals;
   std::map<MVertex*,SVector3> grads,dgrads[3];
 
@@ -26,13 +30,11 @@ class meshMetric: public Field {
   std::map<MVertex*,double> _nodalSizes, _detMetric;
   std::map<MVertex*,fullMatrix<double> > _hessian;
  public:
-  meshMetric(std::vector<MElement*> &, 
-	     simpleFunction<double> *fct, 
-	     meshMetric::MetricComputationTechnique technique, 
-	     double lcmin, double lcmax,
-	     double E = .1, double eps = 1.e-3 );
+  meshMetric(GModel *gm, int technique, 
+	     simpleFunction<double> *fct, std::vector<double> parameters);
+  ~meshMetric();
   inline SMetric3 metricAtVertex (MVertex* v) {return _nodalMetrics[v];}
-  void computeMetric(std::vector<MElement*> &);
+  void computeMetric();
   void computeValues( v2t_cont adj);
   void computeHessian( v2t_cont adj);
   void setAsBackgroundMesh (GModel *gm);
diff --git a/contrib/mmg3d/build/sources/mmg3d1.c b/contrib/mmg3d/build/sources/mmg3d1.c
index a71424f863a933c2d0d90286a82d6418b31e0b35..68e18083a5988d80445dbabb343db99d2db8279f 100644
--- a/contrib/mmg3d/build/sources/mmg3d1.c
+++ b/contrib/mmg3d/build/sources/mmg3d1.c
@@ -155,14 +155,14 @@ MMG_nvoltot=0;
       if ( dd < 5 || dd < 0.05*nd )   break;
       else if ( it > 12 && nd >= na )  break;
     }
-    if ( na + nd > 0 && abs(mesh->info.imprim) > 2 )
+    if ( na + nd > 0 && mesh->info.imprim )
       fprintf(stdout,"     %8d INSERTED   %8d REMOVED   %8d FILTERED\n",
               na,nd,nf);
     }
     while ( na+nd > 0 && ++it < maxtou );
 
 
-  if ( nna+nnd && abs(mesh->info.imprim) < 3 ) {
+  if ( nna+nnd && mesh->info.imprim ) {
     fprintf(stdout,"     %7d INSERTED  %7d REMOVED  %7d FILTERED\n",nna,nnd,nf);
   }
 
diff --git a/contrib/mmg3d/build/sources/mmg3dlib.c b/contrib/mmg3d/build/sources/mmg3dlib.c
index 724671a0d55d65544bd0540c861dfc56896729fd..f613e1e39efa0fa8aef80f143b3687f27d074a83 100644
--- a/contrib/mmg3d/build/sources/mmg3dlib.c
+++ b/contrib/mmg3d/build/sources/mmg3dlib.c
@@ -257,7 +257,7 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
   int k,iadr,i,jj,kk,ii,im;
   double	lambda[3],v[3][3],*mold,*m;
   //fprintf(stdout,"  %s \n", M_STR);
-  fprintf(stdout,"  -- START MMG3d (%d ELEMS)\n", mesh->ne) ;
+  if (opt[6] < 0) fprintf(stdout,"  -- START MMG3d (%d ELEMS)\n", mesh->ne) ;
   if (opt[6] < 0) fprintf(stdout,"  -- MMG3d, Release %s (%s) \n",M_VER,M_REL);
   if (opt[6] < -10) fprintf(stdout,"     Copyright (c) LJLL/IMB, 2010\n");
   if (opt[6] < -10) fprintf(stdout,"    %s\n",COMPIL);
@@ -469,7 +469,7 @@ int MMG_mmg3dlib(int opt[9],MMG_pMesh mesh,MMG_pSol sol) {
     MMG_prilen(mesh,sol);
     MMG_ratio(mesh,sol,NULL);
   }
-  fprintf(stdout,"  -- END MMG3D (%d ELEMS)\n", mesh->ne);
+  if ( MMG_imprim ) fprintf(stdout,"  -- END MMG3D (%d ELEMS)\n", mesh->ne);
   if ( alert )
     fprintf(stdout,"\n  ## WARNING: INCOMPLETE MESH  %d , %d\n",
             mesh->np,mesh->ne);
diff --git a/contrib/mmg3d/build/sources/ratio.c b/contrib/mmg3d/build/sources/ratio.c
index 42e2a0dec11b972e6ae30a389cd5d2a2d62cdc6e..3bfda09cb518732397e111c297c1968c2ce81cd3 100644
--- a/contrib/mmg3d/build/sources/ratio.c
+++ b/contrib/mmg3d/build/sources/ratio.c
@@ -353,7 +353,7 @@ int MMG_ratio(pMesh mesh, pSol sol,char* firaoame) {
   if(inm) fclose(inm);
   
   /* print histo ratio obtained*/
-  if (mesh->info.imprim == 0){
+  if (mesh->info.imprim < 0){
     fprintf(stdout,"        ANISOTROPIC RATIO (MEAN = %6.2g, MAX = %6.2g, MIN = %6.2f)\n",rapavg / rapnum, rapmax, rapmin);
   }
   else if (mesh->info.imprim < 0){
diff --git a/gmshpy/gmshMesh.i b/gmshpy/gmshMesh.i
index 569db4464f0e629651ec1fccd8b0b78965299665..cba7a004aef6c80daebb0d25a25d42941f6bcc13 100644
--- a/gmshpy/gmshMesh.i
+++ b/gmshpy/gmshMesh.i
@@ -12,6 +12,7 @@
   #include "meshPartitionOptions.h"
   #include "Levy3D.h"
   #include "meshPartition.h"
+  #include "meshMetric.h"
 %}
 
 %include "GmshConfig.h"
@@ -22,3 +23,4 @@
 %include "meshPartitionOptions.h"
 %include "Levy3D.h"
 %include "meshPartition.h"
+%include "meshMetric.h"