diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 376f39c7c8584f5f3f2352025e26c759b868d953..b43587128bb71d524859cbfd89d3bb0e1da9300f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
 #include <sstream>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 2aa84adc94e8406dc72bdf1f7957272c4767bc20..fa2d5b4507347c9130b3d039d93fd80339b1f3ab 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
 #include <string.h>
 #include "GmshMessage.h"
 #include "Numeric.h"
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 8546fc8982745a161dc0e18209d983761b530a54..16a29bcb5f278c14171304484a9aea03b11b2718 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
 #include "GmshConfig.h"
 #include "GModel.h"
 #include "MElement.h"
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index dae79acf8c484297d388791e5f0d1570972b949e..c99531621305c5d0051b86cf98412813653c4944 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -39,45 +39,41 @@ void discreteEdge::createTopo()
 
 void discreteEdge::orderMLines()
 {
-  //printf(" *** ORDERING DISCRETE EDGE %d of size %d \n", this->tag(), lines.size());
-
-  std::vector<MLine*> _m ;
+  std::vector<MLine*> _m;
   std::list<MLine*> segments;
 
   //store all lines in a list : segments
   for (unsigned int i = 0; i < lines.size(); i++){
-    //printf("BEFORE ORDERING segment i=%d, vert=%d, %d\n", i, lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
     segments.push_back(lines[i]);
   }
 
   // find a lonly MLine
-  for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+  for (std::list<MLine*>::iterator it = segments.begin(); it != segments.end(); ++it){
     MVertex *vL = (*it)->getVertex(0);
     MVertex *vR = (*it)->getVertex(1);
     std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
-    if (it1==boundv.end())       boundv.insert(std::make_pair(vL,*it));
+    if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
     else  boundv.erase(it1);
     std::map<MVertex*,MLine*>::iterator it2 = boundv.find(vR);
-    if (it2==boundv.end())    boundv.insert(std::make_pair(vR,*it));
-    else    boundv.erase(it2);
+    if (it2 == boundv.end()) boundv.insert(std::make_pair(vR,*it));
+    else boundv.erase(it2);
   }
 
   //find the first MLine and erase it from the list segments
   MLine *firstLine;
-  if (boundv.size()==2){   // non periodic
+  if (boundv.size() == 2){   // non periodic
     firstLine = (boundv.begin())->second;
-    for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+    for (std::list<MLine*>::iterator it = segments.begin(); it != segments.end(); ++it){
       if (*it == firstLine){
         segments.erase(it);
         break;
       }
     }
   }
-  else if (boundv.size()==0) // periodic
-    {
-      firstLine = *(segments.begin());
-      segments.erase(segments.begin());
-    }
+  else if (boundv.size() == 0){ // periodic
+    firstLine = *(segments.begin());
+    segments.erase(segments.begin());
+  }
   else{
     Msg::Error("EdgeCompound %d is wrong (it has %d end points)",tag(),boundv.size());
   }
@@ -87,14 +83,12 @@ void discreteEdge::orderMLines()
   _orientation.push_back(1);
   MVertex *first = _m[0]->getVertex(0);
   MVertex *last = _m[0]->getVertex(1);
-  //printf("first =%d last =%d \n", first->getNum(), last->getNum());
   while (first != last){
     if (segments.empty())break;
     bool found = false;
-    for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end() ; ++it){
+    for (std::list<MLine*>::iterator it = segments.begin(); it != segments.end(); ++it){
       MLine *e = *it;
       if (e->getVertex(0) == last){
-        //printf("orientation 1: beginV=%d \n", e->getVertex(0)->getNum());
         _m.push_back(e);
         segments.erase(it);
         _orientation.push_back(1);
@@ -103,7 +97,6 @@ void discreteEdge::orderMLines()
         break;
       }
       else if (e->getVertex(1) == last){
-        //printf("orientation 0: endV=%d \n", e->getVertex(1)->getNum());
         _m.push_back(e);
         segments.erase(it);
         _orientation.push_back(0);
@@ -134,29 +127,21 @@ void discreteEdge::orderMLines()
   if (_orientation[0] && lines[0]->getVertex(1) != lines[1]->getVertex(1)
       && lines[0]->getVertex(1) != lines[1]->getVertex(0)){
     printf("coucou here \n");
-    for (unsigned int i = 0; i < lines.size(); i++) _orientation[i] = !_orientation[i] ;
+    for (unsigned int i = 0; i < lines.size(); i++) _orientation[i] = !_orientation[i];
   }
-
-//   for (int i=0;i<lines.size();i++){
-//     printf("AFTER ORDERING segment or=%d, vert=%d, %d\n", (int)_orientation[i], lines[i]->getVertex(0)->getNum(), lines[i]->getVertex(1)->getNum() );
-//   }
-
-
-  return;
-
 }
 
 void discreteEdge::setBoundVertices()
 {
 
   if (boundv.size()==2){
-    //printf("Found a regular open Curve \n");
     std::vector<GVertex*> bound_vertices;
-    for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); iter != boundv.end(); iter++){
+    for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); 
+         iter != boundv.end(); iter++){
       MVertex* vE = (iter)->first;
       bool existVertex  = false;
-      for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
-        //printf("Discrete point =%d %d  bound vE=%d\n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum());
+      for(GModel::viter point = model()->firstVertex(); 
+          point != model()->lastVertex(); point++){
         if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
           bound_vertices.push_back((*point));
           existVertex = true;
@@ -164,70 +149,56 @@ void discreteEdge::setBoundVertices()
         }
       }
       if(!existVertex){
-        //printf(" !!! bound vertex %d does not exist, create one \n", vE->getNum());
         GVertex *gvB = new discreteVertex(model(),vE->getNum());
         bound_vertices.push_back(gvB);
         vE->setEntity(gvB);
         gvB->mesh_vertices.push_back(vE);
         gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
         model()->add(gvB);
-     }
-      std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), mesh_vertices.end(), vE) ;
+      }
+      std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), 
+                                                       mesh_vertices.end(), vE);
       if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
-
     }
-
-    //printf("set bound  vertices %d %d size mesh-vertices =%d \n",bound_vertices[0]->tag(),bound_vertices[1]->tag(), mesh_vertices.size());
     v0 = bound_vertices[0];
     v1 = bound_vertices[1];
   }
-  else if (boundv.size()==0){
-    //printf("Found a closed Curve add vertex \n");
+  else if (boundv.size() == 0){
     GVertex* bound_vertex;
     std::vector<MLine*>::const_iterator it = lines.begin();
     MVertex* vE = (*it)->getVertex(0);
     bool existVertex  = false;
-
     for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
-      //printf("Discrete point =%d %d  bound vE=%d \n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum());
-       if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
-         bound_vertex = (*point);
-         existVertex = true;
-         break;
-       }
-     }
-     if(!existVertex){
-  
-       GVertex *gvB = new discreteVertex(model(),vE->getNum());
-       bound_vertex = gvB;
-       gvB->mesh_vertices.push_back(vE);
-       gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
-       model()->add(gvB);
-     }
-
-     std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), mesh_vertices.end(), vE) ;
-     if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
-     
-     //printf("set bound  vertices %d %d \n",bound_vertex->tag(),bound_vertex->tag());
+      if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
+        bound_vertex = (*point);
+        existVertex = true;
+        break;
+      }
+    }
+    if(!existVertex){
+      GVertex *gvB = new discreteVertex(model(),vE->getNum());
+      bound_vertex = gvB;
+      gvB->mesh_vertices.push_back(vE);
+      gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
+      model()->add(gvB);
+    }
+    std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), 
+                                                     mesh_vertices.end(), vE);
+    if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
     v0 = bound_vertex;
     v1 = bound_vertex;
-
   }
 
   v0->addEdge(this);
   v1->addEdge(this);
-
-  return;
 }
 
-
 /*
     +---------------+--------------+----------- ... ----------+
     _pars[0]=0   _pars[1]=1    _pars[2]=2             _pars[N+1]=N+1
 */
 void discreteEdge::parametrize()
 { 
-
   for (unsigned int i = 0; i < lines.size() + 1; i++){
     _pars.push_back(i);
   }
@@ -242,8 +213,8 @@ void discreteEdge::parametrize()
 
   MVertex *vL = getBeginVertex()->mesh_vertices[0];
   int i = 0;
-  for(i = 0 ; i < (int)lines.size() - 1; i++){
-    MVertex *vR ;
+  for(i = 0; i < (int)lines.size() - 1; i++){
+    MVertex *vR;
     if (_orientation[i] == 1 ) vR = lines[i]->getVertex(1);
     else vR = lines[i]->getVertex(0);
     int param = i+1;
@@ -253,54 +224,48 @@ void discreteEdge::parametrize()
     vNEW->setNum(vR->getNum());
     newLines.push_back(new MLine(vL, vNEW));
     _orientation[i] = 1;
-    //printf("i edge =%d , orientation=%d , mesh vertex =%p (%d) newv=%p (%d)  \n", i, _orientation[i], vL, vL->getNum(), vNEW, vNEW->getNum());
     vL = vNEW;
   }
-   MVertex *vR = getEndVertex()->mesh_vertices[0];
-   newLines.push_back(new MLine(vL, vR));
-   _orientation[i] = 1;
-   //printf("i edge =%d , orientation=%d , mesh vertex =%p (%d) newv=%p (%d)  \n", i, _orientation[i], vL, vL->getNum(), vR, vR->getNum());
-
-   mesh_vertices = newVertices;
-   deleteVertexArrays();
-   lines.clear();
-   lines = newLines;
-
-   for(std::list<GFace*>::iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
-     std::vector<MTriangle*> newTriangles;
-     std::vector<MQuadrangle*> newQuadrangles;
-     for (unsigned int i = 0; i < (*iFace)->getNumMeshElements(); ++i){
-       MElement *e = (*iFace)->getMeshElement(i);
-       int N = e->getNumVertices();
-       //MVertex *v[N];
-       std::vector<MVertex *> v(N);
-       for(int j = 0; j < N; j++){
-         v[j] = e->getVertex(j);
-       }
-        //printf("old triangle v0=%p (%d) v1=%p (%d) v2=%p (%d) \n",v[0], v[0]->getNum() , v[1],v[1]->getNum() ,v[2], v[2]->getNum());
-        for (int j = 0; j < N; j++){
-          std::map<MVertex*, MVertex*>::iterator itmap = old2new.find(v[j]);
-          MVertex *vNEW;
-          if (itmap != old2new.end())  {
-            vNEW = itmap->second;
-            v[j]=vNEW;
-          }
+  MVertex *vR = getEndVertex()->mesh_vertices[0];
+  newLines.push_back(new MLine(vL, vR));
+  _orientation[i] = 1;
+
+  mesh_vertices = newVertices;
+  deleteVertexArrays();
+  lines.clear();
+  lines = newLines;
+
+  for(std::list<GFace*>::iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
+    std::vector<MTriangle*> newTriangles;
+    std::vector<MQuadrangle*> newQuadrangles;
+    for (unsigned int i = 0; i < (*iFace)->getNumMeshElements(); ++i){
+      MElement *e = (*iFace)->getMeshElement(i);
+      int N = e->getNumVertices();
+      std::vector<MVertex *> v(N);
+      for(int j = 0; j < N; j++){
+        v[j] = e->getVertex(j);
+      }
+      for (int j = 0; j < N; j++){
+        std::map<MVertex*, MVertex*>::iterator itmap = old2new.find(v[j]);
+        MVertex *vNEW;
+        if (itmap != old2new.end())  {
+          vNEW = itmap->second;
+          v[j]=vNEW;
         }
-        //printf(" new triangle v0=%p (%d) v1=%p (%d) v2=%p (%d) \n",v[0], v[0]->getNum() , v[1],v[1]->getNum() ,v[2], v[2]->getNum());
-        if (N == 3) newTriangles.push_back(new  MTriangle(v[0], v[1], v[2]));
-        else if ( N == 4)  newQuadrangles.push_back(new  MQuadrangle(v[0], v[1], v[2], v[3]));
-
       }
-     (*iFace)->deleteVertexArrays();
-     (*iFace)->triangles.clear();
-     (*iFace)->triangles = newTriangles;
-     (*iFace)->quadrangles.clear();
-     (*iFace)->quadrangles = newQuadrangles;
-   }
-
-
-   //for(std::list<GRegion*>::iterator iRegion = l_regions.begin(); iRegion != l_regions.end(); ++iRegion){
-   for(GModel::riter iRegion = model()->firstRegion(); iRegion != model()->lastRegion(); iRegion++){
+      if (N == 3) newTriangles.push_back(new  MTriangle(v[0], v[1], v[2]));
+      else if ( N == 4)  newQuadrangles.push_back(new  MQuadrangle(v[0], v[1], v[2], v[3]));
+      
+    }
+    (*iFace)->deleteVertexArrays();
+    (*iFace)->triangles.clear();
+    (*iFace)->triangles = newTriangles;
+    (*iFace)->quadrangles.clear();
+    (*iFace)->quadrangles = newQuadrangles;
+  }
+  
+   for(GModel::riter iRegion = model()->firstRegion(); 
+       iRegion != model()->lastRegion(); iRegion++){
      std::vector<MTetrahedron*> newTetrahedra;
      std::vector<MHexahedron*> newHexahedra;
      std::vector<MPrism*> newPrisms;
@@ -308,24 +273,24 @@ void discreteEdge::parametrize()
      for (unsigned int i = 0; i < (*iRegion)->getNumMeshElements(); ++i){
        MElement *e = (*iRegion)->getMeshElement(i);
        int N = e->getNumVertices();
-       //MVertex *v[N];
        std::vector<MVertex *> v(N);
        for(int j = 0; j < N; j++){
          v[j] = e->getVertex(j);
        }
-        for (int j = 0; j < N; j++){
-          std::map<MVertex*, MVertex*>::iterator itmap = old2new.find(v[j]);
-          MVertex *vNEW;
-          if (itmap != old2new.end())  {
-            vNEW = itmap->second;
-            v[j]=vNEW;
-          }
-        }
-        if (N == 4) newTetrahedra.push_back(new  MTetrahedron(v[0], v[1], v[2], v[3]));
-        else if ( N == 5)  newPyramids.push_back(new  MPyramid(v[0], v[1], v[2], v[3], v[4]));
-        else if ( N == 6)  newPrisms.push_back(new  MPrism(v[0], v[1], v[2], v[3], v[4], v[5]));
-        else if ( N == 8)  newHexahedra.push_back(new  MHexahedron(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]));
-      }
+       for (int j = 0; j < N; j++){
+         std::map<MVertex*, MVertex*>::iterator itmap = old2new.find(v[j]);
+         MVertex *vNEW;
+         if (itmap != old2new.end())  {
+           vNEW = itmap->second;
+           v[j]=vNEW;
+         }
+       }
+       if (N == 4) newTetrahedra.push_back(new MTetrahedron(v[0], v[1], v[2], v[3]));
+       else if (N == 5) newPyramids.push_back(new MPyramid(v[0], v[1], v[2], v[3], v[4]));
+       else if (N == 6) newPrisms.push_back(new MPrism(v[0], v[1], v[2], v[3], v[4], v[5]));
+       else if (N == 8) newHexahedra.push_back(new MHexahedron(v[0], v[1], v[2], v[3],
+                                                               v[4], v[5], v[6], v[7]));
+     }
      (*iRegion)->deleteVertexArrays();
      (*iRegion)->tetrahedra.clear();
      (*iRegion)->tetrahedra = newTetrahedra;
@@ -335,28 +300,18 @@ void discreteEdge::parametrize()
      (*iRegion)->prisms = newPrisms;
      (*iRegion)->hexahedra.clear();
      (*iRegion)->hexahedra = newHexahedra;
-
    }
 
    computeNormals();
-
 }
 
 void discreteEdge::computeNormals () const
 {
-  
-//   _normals.clear();
-//   for(std::list<GFace*>::const_iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
-//     GFaceCompound* myFace;
-//     myFace = (GFaceCompound*) *iFace;
-//     myFace->computeNormals(_normals);
-//     printf("coucou \n");
-//   }
-
   _normals.clear();
   double J[3][3];
 
-  for(std::list<GFace*>::const_iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
+  for(std::list<GFace*>::const_iterator iFace = l_faces.begin(); 
+      iFace != l_faces.end(); ++iFace){
     for (unsigned int i = 0; i < (*iFace)->triangles.size(); ++i){
       MTriangle *t = (*iFace)->triangles[i];
       t->getJacobian(0, 0, 0, J);
@@ -372,12 +327,9 @@ void discreteEdge::computeNormals () const
     }
   }
   std::map<MVertex*,SVector3>::iterator itn = _normals.begin();
-  for ( ; itn != _normals.end() ; ++itn){
+  for ( ; itn != _normals.end(); ++itn){
     itn->second.normalize();
-    //printf("NUM = %d xx = %g %g %g  \n", itn->first->getNum(), itn->first->x(), itn->first->y(), itn->first->z() ) ;
-    //printf("normal =%g %g %g \n", itn->second.x(), itn->second.y(), itn->second.z());
   }
-
 }
 
 void discreteEdge::getLocalParameter(const double &t, int &iLine,
@@ -387,8 +339,7 @@ void discreteEdge::getLocalParameter(const double &t, int &iLine,
     double tmin = _pars[iLine];
     double tmax = _pars[iLine+1];
     if (t >= tmin && t <= tmax){
-      tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin)  :  1 - (t-tmin)/(tmax-tmin)  ;
-      //printf("getlocal param t=%g, iLine=%d, tLoc=%g \n", t, iLine, tLoc);
+      tLoc = _orientation[iLine] ? (t-tmin)/(tmax-tmin) : 1 - (t-tmin)/(tmax-tmin);
       return;
     }
   }
@@ -396,7 +347,6 @@ void discreteEdge::getLocalParameter(const double &t, int &iLine,
 
 GPoint discreteEdge::point(double par) const
 {
-
   double tLoc;
   int iEdge;
   getLocalParameter(par,iEdge,tLoc);
@@ -405,23 +355,17 @@ GPoint discreteEdge::point(double par) const
   MVertex *vB = lines[iEdge]->getVertex(0);
   MVertex *vE = lines[iEdge]->getVertex(1);
 
- const bool LINEARMESH = false;
-
+  const bool LINEARMESH = false;
+  
   if (LINEARMESH){
-
     //linear Lagrange mesh
-    //-------------------------
     x = vB->x() + tLoc * (vE->x()- vB->x());
     y = vB->y() + tLoc * (vE->y()- vB->y());
     z = vB->z() + tLoc * (vE->z()- vB->z());
-    
     return GPoint(x,y,z);
   }
   else{
-
     //curved PN triangle
-    //-------------------------
-
     const SVector3 n1 = _normals[vB];
     const SVector3 n2 = _normals[vE];
     
@@ -445,7 +389,6 @@ GPoint discreteEdge::point(double par) const
 
     SPoint3 PP(point.x(), point.y(), point.z());
     return GPoint(PP.x(),PP.y(),PP.z());
-
   }
 
 }
@@ -454,25 +397,22 @@ SVector3 discreteEdge::firstDer(double par) const
 {
   double tLoc;
   int iEdge;
-  getLocalParameter(par,iEdge,tLoc);
+  getLocalParameter(par, iEdge, tLoc);
 
   MVertex *vB = lines[iEdge]->getVertex(0);
   MVertex *vE = lines[iEdge]->getVertex(1);
 
-  double dx,dy,dz;
-  //double dt = sqrt((vE->x()-vB->x())*(vE->x()-vB->x()) + (vE->y()-vB->y())*(vE->y()-vB->y()) + (vE->z()-vB->z())*(vE->z()-vB->z()));
-  //Length of each line of discrete Edge=1.
-  double dt= 1.0;
-  dx = (vE->x() - vB->x() ) /dt;
-  dy = (vE->y() - vB->y() ) /dt;
-  dz = (vE->z() - vB->z() ) /dt;
+  double dx, dy, dz;
+  double dt = 1.0;
+  dx = (vE->x() - vB->x()) / dt;
+  dy = (vE->y() - vB->y()) / dt;
+  dz = (vE->z() - vB->z()) / dt;
 
-  SVector3 der(dx,dy,dz);
+  SVector3 der(dx, dy, dz);
   return der;
 }
 
 Range<double> discreteEdge::parBounds(int i) const
 {
-  //  printf("par bound discreteEdge =%g \n", lines.size());
   return Range<double>(0, lines.size());
 }
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 48f8d67305a3a58e47406d047b79c89d43d0bac6..be66fe48f4cfd57eeb0b47a62d3f2cb180511935 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "discreteFace.h"
@@ -17,11 +18,9 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
   meshStatistics.status = GFace::DONE;    
 }
 
-
 void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges)
 {
-  //find the boundary edges
-
+  // find the boundary edges
   std::list<MEdge> bound_edges;
   for (unsigned int iFace = 0; iFace  < getNumMeshElements() ; iFace++) {
     std::vector<MVertex*> fv;
@@ -29,46 +28,40 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_
     e->getFaceVertices(0, fv);
     for (int iEdge = 0; iEdge < e->getNumEdges(); iEdge++) {
       MEdge tmp_edge =  e->getEdge(iEdge);
-      if (std::find(bound_edges.begin(),bound_edges.end(),tmp_edge) == bound_edges.end())       
+      if (std::find(bound_edges.begin(), bound_edges.end(), tmp_edge) == 
+          bound_edges.end())
         bound_edges.push_back(tmp_edge);
       else
-        bound_edges.erase(std::find(bound_edges.begin(),bound_edges.end(),tmp_edge));
+        bound_edges.erase(std::find(bound_edges.begin(), bound_edges.end(), tmp_edge));
     }
   }
  
-  //for the boundary edges, associate the tag of the current discrete face
-  for (std::list<MEdge>::iterator itv = bound_edges.begin() ; itv != bound_edges.end() ; ++itv){
-    std::map<MEdge, std::vector<int> , Less_Edge >::iterator itmap = map_edges.find(*itv);
-    if (itmap == map_edges.end())   {
+  // for the boundary edges, associate the tag of the current discrete face
+  for (std::list<MEdge>::iterator itv = bound_edges.begin(); 
+       itv != bound_edges.end(); ++itv){
+    std::map<MEdge, std::vector<int>, Less_Edge >::iterator itmap = map_edges.find(*itv);
+    if (itmap == map_edges.end()){
       std::vector<int> tagFaces; 
-      tagFaces.push_back(this->tag());
-      map_edges.insert(std::make_pair(*itv,tagFaces));   
+      tagFaces.push_back(tag());
+      map_edges.insert(std::make_pair(*itv, tagFaces));   
     }
     else{
       std::vector<int> tagFaces = itmap->second;
-      tagFaces.push_back(this->tag());
+      tagFaces.push_back(tag());
       itmap->second = tagFaces;
-      //printf("addface %d %d\n", tagFaces[0], tagFaces[1]);
     }
  }
-
-  //printf( "There are  %d bound msh edges \n ",  map_edges.size());
 }
 
 void discreteFace::setBoundEdges(std::vector<int> tagEdges)
 {
-
-  // printf("***** In discrete Face:  \n");
-
- for (std::vector<int>::iterator it = tagEdges.begin(); it != tagEdges.end(); it++) {
+ for (std::vector<int>::iterator it = tagEdges.begin(); 
+      it != tagEdges.end(); it++){
    GEdge *ge = GModel::current()->getEdgeByTag(abs(*it));
    l_edges.push_back(ge);
    l_dirs.push_back(1);
    ge->addFace(this);
-   //printf("for Face %d add bound edge %d\n",this->tag(), ge->tag() );
  }
-
- //  printf("bound edges =%d \n", edges().size());
 }
 
 GPoint discreteFace::point(double par1, double par2) const 
@@ -79,8 +72,7 @@ GPoint discreteFace::point(double par1, double par2) const
 
 SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
 {
-
-  //Msg::Error("Cannot compute parametric coordinates on discrete face");
+  Msg::Error("Cannot compute parametric coordinates on discrete face");
   return SPoint2();
 }
 
diff --git a/Geo/discreteRegion.cpp b/Geo/discreteRegion.cpp
index 91bddb96a2129e2c8b5276eede1e2cedf1cc16e7..198c41aad772391fe08d6f48c47c1b24165c349b 100644
--- a/Geo/discreteRegion.cpp
+++ b/Geo/discreteRegion.cpp
@@ -21,6 +21,5 @@ void discreteRegion::setBoundFaces()
   for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
     l_faces.push_back(*face);
     (*face)->addRegion(this);
-    printf("face %d \n", (*face)->tag());
   }
 }
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 8846c01aea3512d40cae4ad6e33a24f8fb7b71b3..e8ae1a40340154475bbf3e3a633a4d59d547d518 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -7,6 +7,7 @@
 //   Jonathan Lambrechts
 //
 
+#include <stdlib.h>
 #include <list>
 #include <math.h>
 #include <fstream>
diff --git a/Numeric/gmshLinearSystemCSR.cpp b/Numeric/gmshLinearSystemCSR.cpp
index c9631706c5f4b4672a9361bf9420e2150170c1ad..0e364b437abef6b54391e6e2e468fa56e59a9af5 100644
--- a/Numeric/gmshLinearSystemCSR.cpp
+++ b/Numeric/gmshLinearSystemCSR.cpp
@@ -3,11 +3,12 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "gmshLinearSystemCSR.h"
-#include <stdlib.h>
-#include <string.h>
 
 #define SWAP(a,b)  temp=(a);(a)=(b);(b)=temp;
 #define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;