diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 449b8581225b2816a3537940ba5a07e7a5b8721c..f735595192be05927ba3cab6e15fb8379b6b25ae 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -24,7 +24,6 @@
 #include "discreteEdge.h"
 #include "discreteVertex.h"
 #include "gmshSurface.h"
-#include "Octree.h"
 #include "SmoothData.h"
 #include "Context.h"
 #include "OS.h"
@@ -247,7 +246,7 @@ void GModel::destroyMeshCaches()
   _elementVectorCache.clear();
   _elementMapCache.clear();
   _elementIndexCache.clear();
-  Octree_Delete(_octree);
+  delete _octree;
   _octree = 0;
 }
 
@@ -672,10 +671,9 @@ MElement *GModel::getMeshElementByCoord(SPoint3 &p)
 {
   if(!_octree){
     Msg::Debug("Rebuilding mesh element octree");
-    _octree = buildMElementOctree(this);
+    _octree = new MElementOctree(this);
   }
-  double P[3] = {p.x(), p.y(), p.z()};
-  return (MElement*)Octree_Search(P, _octree);
+  return _octree->find(p.x(), p.y(), p.z());
 }
 
 MVertex *GModel::getMeshVertexByTag(int n)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 19fffce6ac15dea38a428d48a2b9ff494cc9202a..2e2ac3269f90f444e8a9055d8ddc3d6652331a83 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -18,7 +18,6 @@
 #include "SPoint3.h"
 #include "SBoundingBox3d.h"
 
-class Octree;
 class FM_Internals;
 class GEO_Internals;
 class OCC_Internals;
@@ -29,6 +28,7 @@ class CGNSOptions;
 class gLevelset;
 class discreteFace;
 class binding;
+class MElementOctree;
 class GModelFactory;
 
 // A geometric model. The model is a "not yet" non-manifold B-Rep.
@@ -58,7 +58,7 @@ class GModel
   std::multimap<MElement*, short> _ghostCells;
 
   // an octree for fast mesh element lookup
-  Octree *_octree;
+  MElementOctree *_octree;
 
   // Geo (Gmsh native) model internal data
   GEO_Internals *_geo_internals;
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index a1e151adb58d8ac9001b404ffdf15d70b36c77c3..e34fa693fd1c295037dcf6b881ce2281f0f194e3 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -3,10 +3,10 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
-#include "MElement.h"
 #include "GModel.h"
+#include "MElement.h"
+#include "MElementOctree.h"
 #include "Octree.h"
-#include "MVertex.h"
 
 static void MElementBB(void *a, double *min, double *max)
 {
@@ -47,7 +47,7 @@ static int MElementInEle(void *a, double *x)
   return e->isInside(uvw[0], uvw[1], uvw[2]) ? 1 : 0;
 }
 
-Octree *buildMElementOctree(GModel *m)
+MElementOctree::MElementOctree(GModel *m)
 {
   SBoundingBox3d bb = m->bounds();
   double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
@@ -55,18 +55,17 @@ Octree *buildMElementOctree(GModel *m)
                     bb.max().y() - bb.min().y(),
                     bb.max().z() - bb.min().z()};
   const int maxElePerBucket = 100; // memory vs. speed trade-off
-  Octree *_octree = Octree_Create(maxElePerBucket, min, size,
-                                  MElementBB, MElementCentroid, MElementInEle);
+  _octree = Octree_Create(maxElePerBucket, min, size,
+                          MElementBB, MElementCentroid, MElementInEle);
   std::vector<GEntity*> entities;
   m->getEntities(entities);
   for(unsigned int i = 0; i < entities.size(); i++)
     for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
       Octree_Insert(entities[i]->getMeshElement(j), _octree);
   Octree_Arrange(_octree);
-  return _octree;
 }
 
-Octree *buildMElementOctree(std::vector<MElement*> &v)
+MElementOctree::MElementOctree(std::vector<MElement*> &v)
 {
   SBoundingBox3d bb;
   for (unsigned int i = 0; i < v.size(); i++){
@@ -81,10 +80,21 @@ Octree *buildMElementOctree(std::vector<MElement*> &v)
                     bb.max().y() - bb.min().y(),
                     bb.max().z() - bb.min().z()};
   const int maxElePerBucket = 100; // memory vs. speed trade-off
-  Octree *_octree = Octree_Create(maxElePerBucket, min, size,
-                                  MElementBB, MElementCentroid, MElementInEle);
+  _octree = Octree_Create(maxElePerBucket, min, size,
+                          MElementBB, MElementCentroid, MElementInEle);
   for (unsigned int i = 0; i < v.size(); i++)
     Octree_Insert(v[i], _octree);
   Octree_Arrange(_octree);
-  return _octree;
 }
+
+MElementOctree::~MElementOctree()
+{
+  Octree_Delete(_octree);
+}
+
+MElement *MElementOctree::find(double x, double y, double z)
+{
+  double P[3] = {x, y, z};
+  return (MElement*)Octree_Search(P, _octree);
+}
+
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 35606fcd53e72b7da3385ba9f3d1504fb93220e7..262fef6fa3933743deb1fba4b9bad7abefb9cbee 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -12,7 +12,15 @@ class Octree;
 class GModel;
 class MElement;
 
-Octree *buildMElementOctree(GModel *);
-Octree *buildMElementOctree(std::vector<MElement*> &);
+class MElementOctree{
+ private:
+  Octree *_octree;
+ public:
+  MElementOctree(GModel *);
+  MElementOctree(std::vector<MElement*> &);
+  ~MElementOctree();
+  MElement *find(double x, double y, double z);
+  Octree *getInternalOctree(){ return _octree; }
+};
 
 #endif
diff --git a/Geo/MVertexOctree.h b/Geo/MVertexOctree.h
index eb4b5d53c067285f734cc47113fe0c2dbf182a40..79a7806eb6a1aaa6c17cf11652fbba9c2bd81c17 100644
--- a/Geo/MVertexOctree.h
+++ b/Geo/MVertexOctree.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
 #ifndef _MVERTEX_OCTREE_
 #define _MVERTEX_OCTREE_
 
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
index e2c78844a4afe2e3a91a035f99974d8199233105..db90bf13a4f417812af6db326c252cd98d5fabeb 100644
--- a/Geo/STensor3.h
+++ b/Geo/STensor3.h
@@ -30,7 +30,7 @@ class SMetric3 {
       }
     }
   }
-  void setMat (const fullMatrix<double> & mat)
+  void setMat(const fullMatrix<double> & mat)
   {
     for (int i = 0; i < 3; i++)
       for (int j = i; j < 3; j++)
@@ -52,37 +52,37 @@ class SMetric3 {
            const SVector3 &t1,
            const SVector3 &t2,
            const SVector3 &t3)
-    {
-      // M = e^1 * diag * e^1^t
-      // where the elements of diag are l_i = h_i^-2
-      // and the rows of e are the UNIT and ORTHOGONAL directions
-
-      fullMatrix<double> e(3,3);
-      e(0,0) = t1(0); e(0,1) = t1(1); e(0,2) = t1(2);
-      e(1,0) = t2(0); e(1,1) = t2(1); e(1,2) = t2(2);
-      e(2,0) = t3(0); e(2,1) = t3(1); e(2,2) = t3(2);
-      e.transposeInPlace();
-      //      e.invertInPlace();
+  {
+    // M = e^1 * diag * e^1^t
+    // where the elements of diag are l_i = h_i^-2
+    // and the rows of e are the UNIT and ORTHOGONAL directions
     
-      fullMatrix<double> tmp(3,3);
-      tmp(0,0) = l1 * e(0,0);
-      tmp(0,1) = l2 * e(0,1);
-      tmp(0,2) = l3 * e(0,2);
-      tmp(1,0) = l1 * e(1,0);
-      tmp(1,1) = l2 * e(1,1);
-      tmp(1,2) = l3 * e(1,2);
-      tmp(2,0) = l1 * e(2,0);
-      tmp(2,1) = l2 * e(2,1);
-      tmp(2,2) = l3 * e(2,2);
-      
-      e.transposeInPlace();
-
-      _val[0] = tmp(0,0) * e(0,0) + tmp(0,1) * e(1,0) + tmp(0,2) * e(2,0);
-      _val[1] = tmp(1,0) * e(0,0) + tmp(1,1) * e(1,0) + tmp(1,2) * e(2,0);
-      _val[2] = tmp(1,0) * e(0,1) + tmp(1,1) * e(1,1) + tmp(1,2) * e(2,1);
-      _val[3] = tmp(2,0) * e(0,0) + tmp(2,1) * e(1,0) + tmp(2,2) * e(2,0);
-      _val[4] = tmp(2,0) * e(0,1) + tmp(2,1) * e(1,1) + tmp(2,2) * e(2,1);
-      _val[5] = tmp(2,0) * e(0,2) + tmp(2,1) * e(1,2) + tmp(2,2) * e(2,2);
+    fullMatrix<double> e(3,3);
+    e(0,0) = t1(0); e(0,1) = t1(1); e(0,2) = t1(2);
+    e(1,0) = t2(0); e(1,1) = t2(1); e(1,2) = t2(2);
+    e(2,0) = t3(0); e(2,1) = t3(1); e(2,2) = t3(2);
+    e.transposeInPlace();
+    //      e.invertInPlace();
+    
+    fullMatrix<double> tmp(3,3);
+    tmp(0,0) = l1 * e(0,0);
+    tmp(0,1) = l2 * e(0,1);
+    tmp(0,2) = l3 * e(0,2);
+    tmp(1,0) = l1 * e(1,0);
+    tmp(1,1) = l2 * e(1,1);
+    tmp(1,2) = l3 * e(1,2);
+    tmp(2,0) = l1 * e(2,0);
+    tmp(2,1) = l2 * e(2,1);
+    tmp(2,2) = l3 * e(2,2);
+    
+    e.transposeInPlace();
+    
+    _val[0] = tmp(0,0) * e(0,0) + tmp(0,1) * e(1,0) + tmp(0,2) * e(2,0);
+    _val[1] = tmp(1,0) * e(0,0) + tmp(1,1) * e(1,0) + tmp(1,2) * e(2,0);
+    _val[2] = tmp(1,0) * e(0,1) + tmp(1,1) * e(1,1) + tmp(1,2) * e(2,1);
+    _val[3] = tmp(2,0) * e(0,0) + tmp(2,1) * e(1,0) + tmp(2,2) * e(2,0);
+    _val[4] = tmp(2,0) * e(0,1) + tmp(2,1) * e(1,1) + tmp(2,2) * e(2,1);
+    _val[5] = tmp(2,0) * e(0,2) + tmp(2,1) * e(1,2) + tmp(2,2) * e(2,2);
   }
   inline double &operator()(int i, int j)
   {
@@ -175,7 +175,6 @@ SMetric3 interpolation (const SMetric3 &m1,
                         const double v,
                         const double w);
 
-
 // concrete class for general 3x3 matrix
 
 class STensor3 {
@@ -213,7 +212,6 @@ class STensor3 {
     _val[1] = _val[2] = _val[3] = 0.0;
     _val[5] = _val[6] = _val[7] = 0.0;
   }
-
   inline double &operator()(int i, int j)
   {
     return _val[getIndex(i, j)];
@@ -277,14 +275,17 @@ inline double dot(const STensor3 &a, const STensor3 &b)
 }
 
 inline STensor3 operator*(const STensor3 &t, double m)
-{ STensor3 val(t);
-  val*=m;
-  return val; }
+{ 
+  STensor3 val(t);
+  val *= m;
+  return val;
+}
 
 inline STensor3 operator*(double m,const STensor3 &t)
-{ STensor3 val(t);
-  val*=m;
-  return val; }
-
+{ 
+  STensor3 val(t);
+  val *= m;
+  return val;
+}
 
 #endif
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index dda975213897154f2c8271cf22d4344408d9e535..3d36d3bc8d90e83472aba46e6f43bfaef159052f 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -12,12 +12,11 @@
 #include "GFace.h"
 #include "GModel.h"
 #include "Field.h"
-#include "MElementOctree.h"
 #include "MElement.h"
+#include "MElementOctree.h"
 #include "MLine.h"
 #include "MTriangle.h"
 #include "MVertex.h"
-#include "Octree.h"
 
 #if defined(HAVE_SOLVER)
 #include "dofManager.h"
@@ -120,10 +119,7 @@ static SMetric3 metric_based_on_surface_curvature(const GEdge *ge, double u)
       SPoint2 par = ge->reparamOnFace((*it), u, 1);
       SMetric3 m = metric_based_on_surface_curvature (*it, par.x(), par.y());
       if (!count) mesh_size = m;
-      else mesh_size = intersection ( mesh_size, m);    
-      if ( ge->tag() == 197)
-	printf("edge %d face %d g %g %g -> %g %g %g\n",ge->tag(),(*it)->tag(),
-	       m(0,0),m(1,1),m(2,2),mesh_size(0,0),mesh_size(1,1),mesh_size(2,2));
+      else mesh_size = intersection(mesh_size, m);
       count++;
     }
     ++it;
@@ -140,15 +136,17 @@ static SMetric3 metric_based_on_surface_curvature(const GVertex *gv)
     GEdge *_myGEdge = *ite;
     Range<double> bounds = _myGEdge->parBounds(0);
     if (gv == _myGEdge->getBeginVertex())
-      mesh_size = intersection ( mesh_size, 
-				 metric_based_on_surface_curvature (_myGEdge, bounds.low()));    
+      mesh_size = intersection
+        (mesh_size,
+         metric_based_on_surface_curvature(_myGEdge, bounds.low()));
     else
-      mesh_size = intersection ( mesh_size, 
-				 metric_based_on_surface_curvature (_myGEdge, bounds.high()));  }
+      mesh_size = intersection
+        (mesh_size, 
+         metric_based_on_surface_curvature(_myGEdge, bounds.high()));
+  }
   return mesh_size;
 }
 
-
 // the mesh vertex is classified on a model vertex.  we compute the
 // maximum of the curvature of model faces surrounding this point if
 // it is classified on a model edge, we do the same for all model
@@ -184,31 +182,15 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
   return lc;
 }
 
-
 static SMetric3 LC_MVertex_CURV_ANISO(GEntity *ge, double U, double V)
 {
   switch(ge->dim()){
-  case 0:        
-    {
-      SMetric3 m = metric_based_on_surface_curvature((const GVertex *)ge);
-      //      printf("0d\n");
-      return m;
-    }
-    break;
-  case 1:
-    {
-      SMetric3 m = metric_based_on_surface_curvature((const GEdge *)ge,U);
-      //      printf("1d %g %g %g\n",m(0,0),m(1,1),m(2,2));
-      return m;
-    }
-    break;
-  case 2:
-    {
-      return metric_based_on_surface_curvature((const GFace *)ge,U,V);
-    }
-    break;
+  case 0: return metric_based_on_surface_curvature((const GVertex *)ge);
+  case 1: return metric_based_on_surface_curvature((const GEdge *)ge, U);
+  case 2: return metric_based_on_surface_curvature((const GFace *)ge, U, V);
   }
-  Msg::Fatal("curvature control impossible to compute for a volume !");
+  Msg::Error("Curvature control impossible to compute for a volume!");
+  return SMetric3();
 }
 
 // compute the mesh size at a given vertex due to prescribed sizes at
@@ -385,7 +367,7 @@ backgroundMesh::backgroundMesh(GFace *_gf)
     }
     _triangles.push_back(new MTriangle(news[0],news[1],news[2]));
   }
-  _octree = buildMElementOctree(_triangles); 
+  _octree = new MElementOctree(_triangles); 
   if (CTX::instance()->mesh.lcFromPoints)
     propagate1dMesh(_gf);
   else {
@@ -585,12 +567,12 @@ double backgroundMesh::operator() (double u, double v, double w) const
   double uv[3] = {u, v, w};
   double uv2[3];
   //  return 1.0;
-  MElement *e = (MElement*)Octree_Search(uv, _octree);
+  MElement *e = _octree->find(u, v, w);
   if (!e) {
-    Msg::Error("cannot find %g %g",u,v);
+    Msg::Error("cannot find %g %g", u, v);
     return 1.0;
   }
-  e->xyz2uvw (uv,uv2);
+  e->xyz2uvw(uv, uv2);
   std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
   std::map<MVertex*,double>::const_iterator itv2 = _sizes.find(e->getVertex(1));
   std::map<MVertex*,double>::const_iterator itv3 = _sizes.find(e->getVertex(2));
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index f3a78985f55a07fa7da08407c775419fe1648af2..d4b818e7878a0bbb8cb90ecf5caa0e5565847ede 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -10,7 +10,7 @@
 #include "simpleFunction.h"
 #include <vector>
 
-class Octree;
+class MElementOctree;
 class GFace;
 class MElement;
 class MVertex;
@@ -18,7 +18,7 @@ class GEntity;
 
 class backgroundMesh : public simpleFunction<double>
 {
-  Octree *_octree;
+  MElementOctree *_octree;
   std::vector<MVertex*> _vertices;
   std::vector<MElement*> _triangles;
   std::map<MVertex*,double> _sizes;  
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 539e0d3f113cc68b5b30da3b5f1aa78819d47b87..e4102b51a6dde1d7f0031047763226541ca2c004 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -31,7 +31,6 @@
 #include "qualityMeasures.h"
 #include "Field.h"
 #include "OS.h"
-#include "Octree.h"
 #include "MElementOctree.h"
 #include "HighOrder.h"
 #include "meshGEdge.h"
@@ -749,11 +748,10 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
   }
 
   if (Msg::GetVerbosity() == 10){
-    Octree *_octree = buildMElementOctree(gf->model());
+    MElementOctree octree(gf->model());
     doc.Voronoi();
     doc.makePosView("voronoi.pos", gf);
-    doc.printMedialAxis(_octree, "skeleton.pos", gf);
-    Octree_Delete(_octree);
+    doc.printMedialAxis(octree.getInternalOctree(), "skeleton.pos", gf);
   }
 
   gf->triangles.clear();
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 03e78f6e07cb82c91184bb24dfbe626177658af0..c60147a311064d761aa50f9f060dd99c318da4dd 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
 #include <iostream>
 #include "meshGFaceBamg.h"
 #include "meshGFaceLloyd.h"
@@ -13,20 +18,18 @@
 #include <map>
 #include "BackgroundMesh.h"
 
-#ifdef HAVE_BAMG
+#if defined(HAVE_BAMG)
+
 long verbosity = 0;
 #include "Mesh2d.hpp"
 #include "Mesh2.h"
 Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bool);
-//#include "bamg-gmsh.hpp"
-
-static void computeMeshMetricsForBamg (GFace *gf,
-				       int numV,
-				       Vertex2 *bamgVertices,  
-				       double *mm11, 
-				       double *mm12, 
-				       double *mm22,
-				       int iter ){
+
+static void computeMeshMetricsForBamg(GFace *gf, int numV,
+                                      Vertex2 *bamgVertices,  
+                                      double *mm11, double *mm12, double *mm22,
+                                      int iter)
+{
   //  char name[245];
   //  sprintf(name,"bgmBamg-%d-%d.pos",gf->tag(),iter);
   //  if (iter < 2){
@@ -35,11 +38,11 @@ static void computeMeshMetricsForBamg (GFace *gf,
   //  }
 
   fullMatrix<double> J(2,3), JT(3,2),M(3,3),R(2,2),W(2,3);
-  for (int i=0;i<numV;++i){
+  for (int i = 0; i < numV; ++i){
     double u = bamgVertices[i][0];
     double v = bamgVertices[i][1];
-    GPoint gp = gf->point(SPoint2(u,v));
-    SMetric3 m = BGM_MeshMetric(gf,u,v,gp.x(),gp.y(),gp.z());
+    GPoint gp = gf->point(SPoint2(u, v));
+    SMetric3 m = BGM_MeshMetric(gf, u, v, gp.x(), gp.y(), gp.z());
     
     // compute the derivatives of the parametrization
     Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(u, v));
@@ -59,15 +62,15 @@ static void computeMeshMetricsForBamg (GFace *gf,
     mm12[i] = M1.a21;
     mm22[i] = M1.a22;    
   }
-
-
 }
 
-static void meshGFaceBamg_ ( GFace *gf , int iter, bool initialMesh){
+static void meshGFaceBamg_(GFace *gf, int iter, bool initialMesh)
+{
   std::set<MVertex*> all;
   std::map<int,MVertex*> recover;
-  for (int i=0;i<gf->triangles.size();i++){
-    for (int j=0;j<3;j++)all.insert(gf->triangles[i]->getVertex(j));
+  for (unsigned int i = 0; i < gf->triangles.size(); i++){
+    for (unsigned int j = 0; j < 3; j++) 
+      all.insert(gf->triangles[i]->getVertex(j));
   }
   Vertex2 *bamgVertices = new Vertex2[all.size()];
   int index = 0;
@@ -96,8 +99,7 @@ static void meshGFaceBamg_ ( GFace *gf , int iter, bool initialMesh){
   }
 
   Triangle2 *bamgTriangles = new Triangle2[gf->triangles.size()];
-  for (int i=0;i<gf->triangles.size();i++){    
-
+  for (unsigned int i = 0; i < gf->triangles.size(); i++){    
     int nodes [3] = {gf->triangles[i]->getVertex(0)->getIndex(),
 		     gf->triangles[i]->getVertex(1)->getIndex(),
 		     gf->triangles[i]->getVertex(2)->getIndex()};
@@ -113,8 +115,7 @@ static void meshGFaceBamg_ ( GFace *gf , int iter, bool initialMesh){
       nodes[0] = nodes[1];
       nodes[1] = temp;
     }
-    
-    bamgTriangles[i].init (bamgVertices,nodes,gf->tag());
+    bamgTriangles[i].init(bamgVertices, nodes, gf->tag());
   }
   std::list<GEdge*> edges = gf->edges();
   int numEdges = 0;
@@ -126,41 +127,40 @@ static void meshGFaceBamg_ ( GFace *gf , int iter, bool initialMesh){
 
   int count = 0;
   for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
-    for (int i=0; i<(*it)->lines.size();++i){
+    for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
       int nodes [2] = {(*it)->lines[i]->getVertex(0)->getIndex(),
 		       (*it)->lines[i]->getVertex(1)->getIndex()};      
-      bamgBoundary[count].init (bamgVertices,nodes,(*it)->tag());
+      bamgBoundary[count].init (bamgVertices, nodes, (*it)->tag());
       bamgBoundary[count++].lab = count;
     }
   }
   Mesh2 bamgMesh ( all.size(), gf->triangles.size(), numEdges,
-		   bamgVertices, bamgTriangles,bamgBoundary);
-  double *mm11 = new double [all.size()];
-  double *mm12 = new double [all.size()];
-  double *mm22 = new double [all.size()];
+		   bamgVertices, bamgTriangles, bamgBoundary);
+  double *mm11 = new double[all.size()];
+  double *mm12 = new double[all.size()];
+  double *mm22 = new double[all.size()];
   double args[256];
-  computeMeshMetricsForBamg (gf,all.size(),bamgVertices,mm11,mm12,mm22,iter);
-  for (int i=0;i<256;i++)args[i] = -1.1e100;
+  computeMeshMetricsForBamg(gf, all.size(), bamgVertices, mm11, mm12, mm22, iter);
+  for (int i = 0; i < 256; i++) args[i] = -1.1e100;
   Mesh2 *refinedBamgMesh = 0;
   try{
-    refinedBamgMesh = Bamg (&bamgMesh, args, mm11,mm12,mm22, initialMesh);
-    Msg::Info ("bamg succeeded %d vertices %d triangles",
-	       refinedBamgMesh->nv,refinedBamgMesh->nt);
+    refinedBamgMesh = Bamg(&bamgMesh, args, mm11, mm12, mm22, initialMesh);
+    Msg::Info("bamg succeeded %d vertices %d triangles",
+              refinedBamgMesh->nv, refinedBamgMesh->nt);
   }
   catch(...){
-    Msg::Error ("bamg failed");
-    return ;
+    Msg::Error("bamg failed");
+    return;
   }
   delete [] mm11;
   delete [] mm12;
   delete [] mm22;
   std::map<int,MVertex*> yetAnother;
-  for (int i=0;i< refinedBamgMesh->nv;i++){
+  for (int i = 0; i < refinedBamgMesh->nv; i++){
     Vertex2 &v = refinedBamgMesh->vertices[i];
-    //    printf("v.lab = %d\n",v.lab);
     if (i >= nbFixedVertices){
-      GPoint gp = gf->point(SPoint2(v[0],v[1]));
-      MFaceVertex *x = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,v[0],v[1]);
+      GPoint gp = gf->point(SPoint2(v[0], v[1]));
+      MFaceVertex *x = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, v[0], v[1]);
       yetAnother[i] = x;
       gf->mesh_vertices.push_back(x);
     }
@@ -169,38 +169,37 @@ static void meshGFaceBamg_ ( GFace *gf , int iter, bool initialMesh){
     }
   }
 
-  for (int i=0;i<gf->triangles.size();i++){    
+  for (unsigned int i = 0; i < gf->triangles.size(); i++){    
     delete gf->triangles[i];
   }
   gf->triangles.clear();
-  for (int i=0;i< refinedBamgMesh->nt;i++){
+  for (int i = 0; i < refinedBamgMesh->nt; i++){
     Triangle2 &t = refinedBamgMesh->triangles[i];
     Vertex2 &v1 = t[0];
     Vertex2 &v2 = t[1];
     Vertex2 &v3 = t[2];
     gf->triangles.push_back(new MTriangle(yetAnother[(*refinedBamgMesh)(v1)],
 					  yetAnother[(*refinedBamgMesh)(v2)],
-					  yetAnother[(*refinedBamgMesh)(v3)]
-					  ));    
+					  yetAnother[(*refinedBamgMesh)(v3)]));    
   }
-
-
-  if (refinedBamgMesh)delete refinedBamgMesh;
+  
+  if (refinedBamgMesh) delete refinedBamgMesh;
 }
-void meshGFaceBamg ( GFace *gf ){
 
+void meshGFaceBamg(GFace *gf)
+{
   int nT = gf->triangles.size();
   //  meshGFaceBamg_ ( gf , 0, true);
-  for (int i=1;i<14;i++){
+  for (int i = 1; i < 14; i++){
     //    char name[245];
     //    sprintf(name,"hop%d.msh",i);
     //    GModel::current()->writeMSH(name);
-    meshGFaceBamg_ ( gf , i, false);
+    meshGFaceBamg_(gf, i, false);
     //    sprintf(name,"hap%d.msh",i);
     //    GModel::current()->writeMSH(name);
     
     int nTnow = gf->triangles.size();
-    if (fabs((double)(nTnow - nT)) < 0.01 * nT)break;
+    if (fabs((double)(nTnow - nT)) < 0.01 * nT) break;
     nT = nTnow;
   }
   //  lloydAlgorithm lloyd (20);
@@ -208,7 +207,11 @@ void meshGFaceBamg ( GFace *gf ){
 }
 
 #else
-void meshGFaceBamg ( GFace *gf ){
-  Msg::Fatal("Gmsh msust be compiled with the Bidimensional Anisotropic Mesh Generator (Bamg) in order to support that option");
+
+void meshGFaceBamg(GFace *gf)
+{
+  Msg::Error("Gmsh msust be compiled with the Bidimensional Anisotropic Mesh "
+             "Generator (Bamg) in order to support that option");
 }
+
 #endif
diff --git a/Mesh/meshGFaceBamg.h b/Mesh/meshGFaceBamg.h
index dc57394256970a4303defc0f0249838f9e319267..0a90874e980304645d53b1823668001acfd471f0 100644
--- a/Mesh/meshGFaceBamg.h
+++ b/Mesh/meshGFaceBamg.h
@@ -1,6 +1,12 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
 #ifndef _MESHGFACE_BAMG_
 #define _MESHGFACE_BAMG_
+
 class GFace;
-void meshGFaceBamg ( GFace *gf );
+void meshGFaceBamg(GFace *gf);
 
 #endif
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 8fbccdf9ff2d4f31957520d6cfb8f586439a9dcc..6e159f8977ab80f39daf64e9b7bec0721b551276 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -530,8 +530,7 @@ void DocRecord::ConvertDListToVoronoiData()
 void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
 {
   if (!_adjacencies){
-    printf("no adjacencies were created\n");
-    throw;
+    Msg::Error("No adjacencies were created");
   }
   const int n = _adjacencies[pt].t_length;
   for(int j = 0; j < n; j++) {