diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 3847d4899c38fe307864b56773ee7317484de3ed..03f050de6c6917c584e3059fa9ef20464e1d83f4 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.98 2008-02-16 20:42:40 geuzaine Exp $
+// $Id: BDS.cpp,v 1.99 2008-02-16 21:25:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -40,14 +40,14 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
   while(tit != tite) {
     BDS_Point *pts[4];
     (*tit)->getNodes(pts);
-    if (param)
+    if(param)
       fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
 	      pts[0]->u, pts[0]->v, 0.0,
 	      pts[1]->u, pts[1]->v, 0.0,
 	      pts[2]->u, pts[2]->v, 0.0,
 	      pts[0]->iD, pts[1]->iD, pts[2]->iD);
     else{
-      if (!gf)
+      if(!gf)
 	fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
 		pts[0]->X, pts[0]->Y, pts[0]->Z,
 		pts[1]->X, pts[1]->Y, pts[1]->Z,
@@ -68,47 +68,47 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
   fclose(f);
 }
 
-BDS_Vector::BDS_Vector(const BDS_Point & p2, const BDS_Point & p1)
+BDS_Vector::BDS_Vector(const BDS_Point &p2, const BDS_Point &p1)
   :x(p2.X - p1.X), y(p2.Y - p1.Y), z(p2.Z - p1.Z)
 {
 }
 
-void vector_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3,
+void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
                      double c[3])
 {
-  double a[3] = { p1->X - p2->X, p1->Y - p2->Y, p1->Z - p2->Z };
-  double b[3] = { p1->X - p3->X, p1->Y - p3->Y, p1->Z - p3->Z };
+  double a[3] = {p1->X - p2->X, p1->Y - p2->Y, p1->Z - p2->Z};
+  double b[3] = {p1->X - p3->X, p1->Y - p3->Y, p1->Z - p3->Z};
   c[2] = a[0] * b[1] - a[1] * b[0];
   c[1] = -a[0] * b[2] + a[2] * b[0];
   c[0] = a[1] * b[2] - a[2] * b[1];
 }
 
-void vector_triangle_parametric(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3,
+void vector_triangle_parametric(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
 				double &c)
 {
-  double a[2] = { p1->u - p2->u, p1->v - p2->v };
-  double b[2] = { p1->u - p3->u, p1->v - p3->v };
+  double a[2] = {p1->u - p2->u, p1->v - p2->v};
+  double b[2] = {p1->u - p3->u, p1->v - p3->v};
   c = a[0] * b[1] - a[1] * b[0];
 }
 
-void normal_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3,
+void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
                      double c[3])
 {
   vector_triangle(p1, p2, p3, c);
   norme(c);
 }
 
-double surface_triangle(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
+double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
 {
   double c[3];
   vector_triangle(p1, p2, p3, c);
   return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
 }
 
-double surface_triangle_param(BDS_Point * p1, BDS_Point * p2, BDS_Point * p3)
+double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
 {
   double c;
-  vector_triangle_parametric (p1, p2, p3, c);
+  vector_triangle_parametric(p1, p2, p3, c);
   return (0.5 * c);
 }
 
@@ -162,14 +162,14 @@ BDS_Point *BDS_Mesh::find_point(int p)
   BDS_Point P(p);
   std::set < BDS_Point *, PointLessThan >::iterator it = points.find(&P);
   if(it != points.end())
-    return (BDS_Point *) (*it);
+    return (BDS_Point*)(*it);
   else
     return 0;
 }
 
 BDS_Edge *BDS_Mesh::find_edge(BDS_Point *p, int num2)
 {
-  std::list < BDS_Edge * >::iterator eit = p->edges.begin();
+  std::list<BDS_Edge*>::iterator eit = p->edges.begin();
   while(eit != p->edges.end()) {
     if((*eit)->p1 == p && (*eit)->p2->iD == num2)
       return (*eit);
@@ -205,7 +205,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
 //   double rq1 = gmsh::orient2d(q1, q2, p1);
 //   double rq2 = gmsh::orient2d(q1, q2, p2);
 
-//   if (rp1*rp2<=0 && rq1*rq2<=0){
+//   if(rp1*rp2<=0 && rq1*rq2<=0){
 // //      printf("p1 %22.15E %22.15E %22.15E %22.15E \n",x1,y1,x2,y2);
 // //      printf("p2 %22.15E %22.15E %22.15E %22.15E \n",x3,y3,x4,y4);
 // //      printf("or %22.15E %22.15E %22.15E %22.15E\n",rp1,rp2,rq1,rq2);
@@ -231,14 +231,14 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
 BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2r, 
 				 std::set<EdgeToRecover> *not_recovered)
 {
-  BDS_Edge *e = find_edge (num1, num2);
+  BDS_Edge *e = find_edge(num1, num2);
 
-  if (e) return e;
+  if(e) return e;
 
   BDS_Point *p1 = find_point(num1);
   BDS_Point *p2 = find_point(num2);
   
-  if (!p1 || !p2) throw;
+  if(!p1 || !p2) throw;
 
   Msg(DEBUG, " edge %d %d has to be recovered", num1, num2);
   
@@ -273,7 +273,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
       ++it;
     }
 
-//   if (ix > 300){
+//   if(ix > 300){
 //     Msg(WARNING," edge %d %d cannot be recovered after %d iterations, trying again",
 // 	  num1, num2, ix);	
 //     ix = 0;
@@ -281,7 +281,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, std::set<EdgeToRecover> *e2
 //   printf("%d %d\n",intersected.size(),ix);
       
     if(!intersected.size() || ix > 1000){
-      BDS_Edge *eee = find_edge (num1, num2);
+      BDS_Edge *eee = find_edge(num1, num2);
       if(!eee){
 	outputScalarField(triangles, "debugp.pos", 1);
 	outputScalarField(triangles, "debugr.pos", 0);
@@ -400,7 +400,7 @@ BDS_Face *BDS_Mesh::add_quadrangle(BDS_Edge *e1, BDS_Edge *e2,
   return t;
 }
 
-void BDS_Mesh::del_face(BDS_Face * t)
+void BDS_Mesh::del_face(BDS_Face *t)
 {
   t->e1->del(t);
   t->e2->del(t);
@@ -409,14 +409,14 @@ void BDS_Mesh::del_face(BDS_Face * t)
   t->deleted = true;
 }
 
-void BDS_Mesh::del_edge(BDS_Edge * e)
+void BDS_Mesh::del_edge(BDS_Edge *e)
 {
   e->p1->del(e);
   e->p2->del(e);
   e->deleted = true;
 }
 
-void BDS_Mesh::del_point(BDS_Point * p)
+void BDS_Mesh::del_point(BDS_Point *p)
 {
   points.erase(p);
   delete p;
@@ -427,7 +427,7 @@ void BDS_Mesh::add_geom(int p1, int p2)
   geom.insert(new BDS_GeomEntity(p1, p2));
 }
 
-void BDS_Edge::oppositeof(BDS_Point * oface[2]) const
+void BDS_Edge::oppositeof(BDS_Point *oface[2]) const
 {
   oface[0] = oface[1] = 0;
   if(faces(0)) {
@@ -457,10 +457,10 @@ BDS_GeomEntity *BDS_Mesh::get_geom(int p1, int p2)
   BDS_GeomEntity ge(p1, p2);
   std::set < BDS_GeomEntity *, GeomLessThan >::iterator it = geom.find(&ge);
   if(it == geom.end()) return 0;
-  return (BDS_GeomEntity *) (*it);
+  return (BDS_GeomEntity*)(*it);
 }
 
-void recur_tag(BDS_Face * t, BDS_GeomEntity * g)
+void recur_tag(BDS_Face *t, BDS_GeomEntity *g)
 {
   if(!t->g) {
     t->g = g;
@@ -564,7 +564,7 @@ bool BDS_Mesh::split_face(BDS_Face *f, BDS_Point *mid)
   return true;
 }
 
-bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
+bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
 {
   /*
         p1
@@ -597,7 +597,6 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
         orientation = 1;
       else
         orientation = -1;
-
       break;
     }
   }
@@ -678,7 +677,7 @@ bool BDS_Mesh::split_edge(BDS_Edge * e, BDS_Point *mid)
 bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2,
 					   BDS_Point *_q1, BDS_Point *_q2) const
 {  
-  if (!testSmallTriangles){
+  if(!testSmallTriangles){
     double p1 [2] = {_p1->u,_p1->v};
     double p2 [2] = {_p2->u,_p2->v};
     double op1[2] = {_q1->u,_q1->v};
@@ -695,8 +694,8 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2,
   double s2 = fabs(surface_triangle_param(_p1, _p2, _q2)); 
   double s3 = fabs(surface_triangle_param(_p1, _q1, _q2)); 
   double s4 = fabs(surface_triangle_param(_p2, _q1, _q2)); 
-  if (fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s1 + s2)) return false;
-  if (s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
+  if(fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s1 + s2)) return false;
+  if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
     return false;
   return true;
 }
@@ -706,7 +705,7 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_P
 					   BDS_Point *_op1, BDS_Point *_op2, BDS_Point *_op3,
 					   BDS_Point *_oq1, BDS_Point *_oq2, BDS_Point *_oq3) const
 {
-  if (!testQuality) return true;
+  if(!testQuality) return true;
   double n[3], q[3], on[3], oq[3];
   normal_triangle(_p1, _p2, _p3, n); 
   normal_triangle(_q1, _q2, _q3, q); 
@@ -725,17 +724,17 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_P
   double mina = std::min(qa1,qa2);
   double minb = std::min(qb1,qb2);
 
-  // if (cosnq < .3 && cosonq > .5 && minb > .1)
+  // if(cosnq < .3 && cosonq > .5 && minb > .1)
   //   printf("mina = %g minb = %g cos %g %g\n",mina,minb,cosnq,cosonq);
   
-  if (cosnq < .3 && cosonq > .5 && minb > .1) return true;
+  if(cosnq < .3 && cosonq > .5 && minb > .1) return true;
   
-  if (minb > mina) return true;
-  //  if (mina > minb && cosnq <= cosonq)return true;
+  if(minb > mina) return true;
+  //  if(mina > minb && cosnq <= cosonq)return true;
   return false;
 }
 
-void swap_config(BDS_Edge * e, 
+void swap_config(BDS_Edge *e, 
 		 BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
 		 BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
 		 BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
@@ -848,21 +847,21 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest)
   }
 
   if(orientation == 1) {
-    if (!theTest(p1, p2, op[0],
-		 p2, p1, op[1],
-		 p1, op[1], op[0],
-		 op[1], p2, op[0]))
+    if(!theTest(p1, p2, op[0],
+		p2, p1, op[1],
+		p1, op[1], op[0],
+		op[1], p2, op[0]))
       return false;
   }
   else{
-    if (!theTest(p2, p1, op[0],
-		 p1, p2, op[1],
-		 p1, op[0], op[1],
-		 op[1], op[0], p2))
+    if(!theTest(p2, p1, op[0],
+		p1, p2, op[1],
+		p1, op[0], op[1],
+		op[1], op[0], p2))
       return false;
   }
   
-  if (!theTest(p1, p2, op[0], op[1]))
+  if(!theTest(p1, p2, op[0], op[1]))
     return false;
 
   BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
@@ -913,15 +912,15 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest)
 int BDS_Edge::numTriangles() const
 {
   int NT = 0;
-  for (unsigned int i = 0; i < _faces.size(); i++) 
-    if (faces(i)->numEdges() == 3) NT++ ;
+  for(unsigned int i = 0; i < _faces.size(); i++) 
+    if(faces(i)->numEdges() == 3) NT++ ;
   return NT;
 }
 
 // This function does actually the swap without taking into account
 // the feasability of the operation. Those conditions have to be
 // taken into account before doing the edge swap
-bool BDS_Mesh::recombine_edge(BDS_Edge * e)
+bool BDS_Mesh::recombine_edge(BDS_Edge *e)
 {
   /*
         p1
@@ -937,7 +936,6 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
    */
 
   // we test if the edge is deleted
-  //  return false;
   if(e->deleted)
     return false;
   if(e->numfaces() != 2 || e->numTriangles() != 2)
@@ -945,24 +943,21 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
   if(e->g && e->g->classif_degree == 1)
     return false;
 
-  BDS_Point *op[2];
   BDS_Point *p1 = e->p1;
   BDS_Point *p2 = e->p2;
+  BDS_Point *op[2];
   e->oppositeof(op);
 
   BDS_Point *pts1[4];
   e->faces(0)->getNodes(pts1);
 
-  //  FIXME  !!!!!!!!!!!!!!!!!
-  //  should ensure that orientation is unchanged
-
   BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
   BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
   BDS_Edge *p1_op2 = find_edge(p1, op[1], e->faces(1));
   BDS_Edge *op2_p2 = find_edge(op[1], p2, e->faces(1));
 
-  BDS_GeomEntity *g=0;
-  if(e->faces(0)) {
+  BDS_GeomEntity *g = 0;
+  if(e->faces(0)){
     g = e->faces(0)->g;
     del_face(e->faces(0));
   }
@@ -972,9 +967,24 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
   }
   del_edge(e);
 
-  BDS_Face *f = add_quadrangle (p1_op1, op1_p2, op2_p2, p1_op2);
-  f->g = g;
+  int orientation = 0;
+  for(int i = 0; i < 3; i++) {
+    if(pts1[i] == p1) {
+      if(pts1[(i + 1) % 3] == p2)
+        orientation = 1;
+      else
+        orientation = -1;
+      break;
+    }
+  }
 
+  BDS_Face *f;
+  if(orientation < 0)
+    f = add_quadrangle(p1_op1, op1_p2, op2_p2, p1_op2);
+  else
+    f = add_quadrangle(p1_op1, p1_op2, op2_p2, op1_p2);
+  f->g = g;
+  
   p1->config_modified = true;
   p2->config_modified = true;
   op[0]->config_modified = true;
@@ -983,7 +993,7 @@ bool BDS_Mesh::recombine_edge(BDS_Edge * e)
   return true;
 }
 
-bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
+bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
 {
   if(e->numfaces() != 2)
     return false;
@@ -997,7 +1007,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
       return false;
   }
 
-  std::list < BDS_Face * >t;
+  std::list<BDS_Face*> t;
   BDS_Point *o = e->othervertex(p);
 
   //  if(o->g != p->g)
@@ -1006,19 +1016,19 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
   // printf("collapsing an edge :");
   // print_edge(e);
 
-  static BDS_Point* pt[3][1024];
+  static BDS_Point *pt[3][1024];
   static BDS_GeomEntity *gs[1024];
   static int ept[2][1024];
   static BDS_GeomEntity *egs[1024]; 
   int nt = 0;
   {
     p->getTriangles(t);
-    std::list < BDS_Face * >::iterator it = t.begin();
-    std::list < BDS_Face * >::iterator ite = t.end();
+    std::list<BDS_Face*>::iterator it = t.begin();
+    std::list<BDS_Face*>::iterator ite = t.end();
     while(it != ite) {
       BDS_Face *t = *it;
       if(t->e1 != e && t->e2 != e && t->e3 != e) {
-	if (!test_move_point_parametric_triangle(p, o->u, o->v, t))
+	if(!test_move_point_parametric_triangle(p, o->u, o->v, t))
 	  return false;
         gs[nt] = t->g;
 	BDS_Point *pts[4];
@@ -1040,8 +1050,8 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge * e, BDS_Point * p)
   }
 
   {
-    std::list < BDS_Face * >::iterator it = t.begin();
-    std::list < BDS_Face * >::iterator ite = t.end();
+    std::list<BDS_Face*>::iterator it = t.begin();
+    std::list<BDS_Face*>::iterator ite = t.end();
 
     while(it != ite) {
       BDS_Face *t = *it;
@@ -1101,19 +1111,19 @@ bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_F
 
   double area_init = fabs(a[0] * b[1] - a[1] * b[0]);
 
-  if (area_init == 0.0) return true;
+  if(area_init == 0.0) return true;
 
   double ori_init = gmsh::orient2d(pa, pb, pc);
 
-  if (p == pts[0]){ 
+  if(p == pts[0]){ 
     pa[0] = u; 
     pa[1] = v; 
   }
-  else if (p == pts[1]){
+  else if(p == pts[1]){
     pb[0] = u;
     pb[1] = v;
   }
-  else if (p == pts[2]){
+  else if(p == pts[2]){
     pc[0] = u;
     pc[1] = v;
   }
@@ -1121,17 +1131,14 @@ bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_F
     return false;
   
   double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
-  if (area_final < 0.1 * area_init) return false;
+  if(area_final < 0.1 * area_init) return false;
   double ori_final = gmsh::orient2d(pa, pb, pc);
   // allow to move a point when a triangle was flat
   return ori_init*ori_final > 0;
 }
 
-/**
-   d^2_i = (x^2_i - x)^T M (x_i - x)  
-         = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2	 
-	~:-)
-*/
+// d^2_i = (x^2_i - x)^T M (x_i - x)  
+//       = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2	 
 
 struct smoothVertexData{
   BDS_Point *p;
@@ -1192,7 +1199,7 @@ double smooth_obj(double U, double V, void *data)
 				      svd->scalu, svd->scalv, svd->gf); 
 }
 
-void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv)
+void optimize_vertex_position(GFace *GF, BDS_Point *data, double su, double sv)
 {
 #ifdef HAVE_GSL
   if(data->g && data->g->classif_degree <= 1) return;
@@ -1205,14 +1212,14 @@ void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv)
   double U = data->u, V = data->v, val;
 
   val = smooth_obj(U, V, &vd);
-  if (val < -.90) return;
+  if(val < -.90) return;
 
   minimize_2(smooth_obj, deriv_smoothing_objective_function, &vd, 5, U,V,val);
   std::list<BDS_Face*>::iterator it = vd.ts.begin();
   std::list<BDS_Face*>::iterator ite = vd.ts.end();
   while(it != ite) {
     BDS_Face *t = *it;
-    if (!test_move_point_parametric_triangle(data, U, V, t)){
+    if(!test_move_point_parametric_triangle(data, U, V, t)){
       return;          
     }
     ++it;
@@ -1229,12 +1236,12 @@ void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv)
 
 bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
 {
-  if (!p->config_modified) return false;
+  if(!p->config_modified) return false;
   if(p->g && p->g->classif_degree <= 1) return false;
 
-  std::list < BDS_Edge * >::iterator eit = p->edges.begin();
+  std::list<BDS_Edge*>::iterator eit = p->edges.begin();
   while(eit != p->edges.end()) {
-    if ((*eit)->numfaces() == 1) return false;
+    if((*eit)->numfaces() == 1) return false;
     eit++;
   }
 
@@ -1290,10 +1297,10 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
     p->v = oldV;
     double sold = fabs(surface_triangle_param(n[0], n[1], n[2])); 
     s2 += sold;
-    //    printf("%22.15E %22.15E\n",snew,sold);
-    if (snew < .1 * sold) return false;
+    // printf("%22.15E %22.15E\n",snew,sold);
+    if(snew < .1 * sold) return false;
 
-    if (1 || test_quality){
+    if(1 || test_quality){
       p->X = gp.x();
       p->Y = gp.y();
       p->Z = gp.z();
@@ -1305,16 +1312,17 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
       p->Z = oldZ;
       normal_triangle(n[0], n[1], n[2], norm2);
       oldWorst = std::min(oldWorst, qmTriangle(*it, QMTRI_RHO));
-      double ps; prosca (norm1, norm2, &ps);
-      if (ps < .5) return false;
+      double ps; 
+      prosca(norm1, norm2, &ps);
+      if(ps < .5) return false;
     }
     ++it;
   }
   
-  //  printf("%22.15E %22.15E %22.15E\n",s1,s2,fabs(s2-s1));
-  if (fabs(s2-s1) > 1.e-14 * (s2 + s1)) return false;
+  // printf("%22.15E %22.15E %22.15E\n",s1,s2,fabs(s2-s1));
+  if(fabs(s2-s1) > 1.e-14 * (s2 + s1)) return false;
   
-  if (test_quality && newWorst < oldWorst){
+  if(test_quality && newWorst < oldWorst){
     return false;
   }
   
@@ -1332,11 +1340,11 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
   return true;
 }
 
-bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
+bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
 {
 
-  if (!p->config_modified)return false;
-  if (p->g && p->g->classif_degree <= 1)
+  if(!p->config_modified)return false;
+  if(p->g && p->g->classif_degree <= 1)
     return false;
   
   double U = 0;
@@ -1344,9 +1352,9 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
   double tot_length = 0; 
   double LC = 0;
 
-  std::list < BDS_Edge * >::iterator eit = p->edges.begin();
+  std::list<BDS_Edge*>::iterator eit = p->edges.begin();
   while(eit != p->edges.end()) {
-    if ((*eit)->numfaces() == 1) return false;
+    if((*eit)->numfaces() == 1) return false;
     BDS_Point *op = ((*eit)->p1 == p) ? (*eit)->p2 : (*eit)->p1;
     const double l_e = (*eit)->length();     
     U += op->u * l_e; 
@@ -1366,7 +1374,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point * p, GFace *gf)
   std::list<BDS_Face*>::iterator ite = ts.end();
   while(it != ite) {
     BDS_Face *t = *it;
-    if (!test_move_point_parametric_triangle(p, U, V, t))
+    if(!test_move_point_parametric_triangle(p, U, V, t))
       return false;
     ++it;
   }
@@ -1409,43 +1417,43 @@ recombine_T::recombine_T(const BDS_Edge *_e)
   BDS_Vector v2(*op[0], *p2);
   BDS_Vector v3(*p2, *op[1]);
   BDS_Vector v4(*op[1], *p1);
-  angle = fabs(90.-v1.angle_deg(v2));
-  angle = std::max(fabs(90.-v2.angle_deg(v3)), angle);
-  angle = std::max(fabs(90.-v3.angle_deg(v4)), angle);
-  angle = std::max(fabs(90.-v4.angle_deg(v1)), angle);
+  angle = fabs(90. - v1.angle_deg(v2));
+  angle = std::max(fabs(90. - v2.angle_deg(v3)), angle);
+  angle = std::max(fabs(90. - v3.angle_deg(v4)), angle);
+  angle = std::max(fabs(90. - v4.angle_deg(v1)), angle);
 }
 
-void BDS_Mesh::recombineIntoQuads (const double angle_limit, GFace *gf)
+void BDS_Mesh::recombineIntoQuads(const double angle_limit, GFace *gf)
 {
   
-  Msg(INFO,"Recombining triangles for surface %d",gf->tag());  
-  for (int i = 0; i < 5 ; i++){
+  Msg(INFO,"Recombining triangles for surface %d", gf->tag());  
+  for(int i = 0; i < 5; i++){
     bool rec = false;
     std::set<recombine_T> _pairs;
     
-    for(std::list < BDS_Edge * >::const_iterator it = edges.begin();
+    for(std::list<BDS_Edge*>::const_iterator it = edges.begin();
 	it != edges.end(); ++it){
-      if (!(*it)->deleted && (*it)->numfaces () == 2)
+      if(!(*it)->deleted && (*it)->numfaces () == 2)
 	_pairs.insert(recombine_T(*it));
     }  
     
     std::set<recombine_T>::iterator itp = _pairs.begin();
     
-    while (itp != _pairs.end()){
-      rec |= recombine_edge ((BDS_Edge*)itp->e);
+    while(itp != _pairs.end()){
+      rec |= recombine_edge((BDS_Edge*)itp->e);
       ++itp;
     }
-    
-    if (!rec) break;
 
-    std::set<BDS_Point*,PointLessThan>::iterator itpt = points.begin();
-    while (itpt != points.end()){
-      smooth_point_parametric(*itpt,gf);
+    if(!rec) break;
+
+    std::set<BDS_Point*, PointLessThan>::iterator itpt = points.begin();
+    while(itpt != points.end()){
+      smooth_point_parametric(*itpt, gf);
       ++itpt;
     }
   }
 }
 
-void FullQuadMesh ()
+void FullQuadMesh()
 {
 }
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index c8882ae66cf363d09dba0110d942e42369994dda..ade3620dc33c30b3ee0e80cfec88914bf53184d8 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -48,12 +48,12 @@ void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
 void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]); 
 double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3); 
 double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3); 
-void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv);
-void swap_config(BDS_Edge * e, 
-		 BDS_Point **p11,BDS_Point **p12,BDS_Point **p13,
-		 BDS_Point **p21,BDS_Point **p22,BDS_Point **p23,
-		 BDS_Point **p31,BDS_Point **p32,BDS_Point **p33,
-		 BDS_Point **p41,BDS_Point **p42,BDS_Point **p43);
+void optimize_vertex_position(GFace *GF, BDS_Point *data, double su, double sv);
+void swap_config(BDS_Edge *e, 
+		 BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
+		 BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
+		 BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
+		 BDS_Point **p41, BDS_Point **p42, BDS_Point **p43);
 
 
 class BDS_GeomEntity
@@ -62,7 +62,7 @@ public:
 
   int classif_tag;
   int classif_degree;
-    inline bool operator < (const BDS_GeomEntity & other) const
+  inline bool operator < (const BDS_GeomEntity & other) const
   {
     if(classif_degree < other.classif_degree)return true;
     if(classif_degree > other.classif_degree)return false;
@@ -93,74 +93,74 @@ public:
     if(z - o.z  > t ) return true;
     return false;
   }
-  BDS_Vector operator + (const  BDS_Vector &v)
+  BDS_Vector operator + (const BDS_Vector &v)
   {
-    return BDS_Vector(x+v.x,y+v.y,z+v.z);
+    return BDS_Vector(x + v.x, y + v.y, z + v.z);
   }
-  BDS_Vector operator - (const  BDS_Vector &v)
+  BDS_Vector operator - (const BDS_Vector &v)
   {
-    return BDS_Vector(x-v.x,y-v.y,z-v.z);
+    return BDS_Vector(x - v.x, y - v.y, z - v.z);
   }
   inline BDS_Vector operator % (const BDS_Vector &other) const
   {
-    return BDS_Vector(y*other.z-z*other.y,
-		      z*other.x-x*other.z,
-		      x*other.y-y*other.x);
+    return BDS_Vector(y * other.z - z * other.y,
+		      z * other.x - x * other.z,
+		      x * other.y - y * other.x);
   }
-  BDS_Vector& operator += (const  BDS_Vector &v)
+  BDS_Vector& operator += (const BDS_Vector &v)
   {
-    x+=v.x;
-    y+=v.y;
-    z+=v.z;
+    x += v.x;
+    y += v.y;
+    z += v.z;
     return *this;
   }
-  BDS_Vector& operator *= (const  double &v)
+  BDS_Vector& operator *= (const double &v)
   {
-    x*=v;
-    y*=v;
-    z*=v;
+    x *= v;
+    y *= v;
+    z *= v;
     return *this;
   }
-  BDS_Vector& operator /= (const  double &v)
+  BDS_Vector& operator /= (const double &v)
   {
     x/=v;
     y/=v;
     z/=v;
     return *this;
   }
-  BDS_Vector operator / (const  double &v)
+  BDS_Vector operator / (const double &v)
   {
-    return BDS_Vector(x/v,y/v,z/v);
+    return BDS_Vector(x / v, y / v, z / v);
   }
-  BDS_Vector operator * (const  double &v)
+  BDS_Vector operator * (const double &v)
   {
-    return BDS_Vector(x*v,y*v,z*v);
+    return BDS_Vector(x * v, y * v, z * v);
   }
-  double angle(const  BDS_Vector &v) const
+  double angle(const BDS_Vector &v) const
   {
-    double a[3] = { x ,  y ,  z };
-    double b[3] = { v.x ,  v.y ,  v.z };
+    double a[3] = {x, y, z};
+    double b[3] = {v.x, v.y, v.z};
     double c[3];
     c[2] = a[0] * b[1] - a[1] * b[0];
     c[1] = -a[0] * b[2] + a[2] * b[0];
     c[0] = a[1] * b[2] - a[2] * b[1];
-    double cosa = a[0]*b[0] +a[1]*b[1] +a[2]*b[2];
-    double sina = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
-    double ag = atan2(sina,cosa);
+    double cosa = a[0] * b[0] + a[1] * b[1] +a[2] * b[2];
+    double sina = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
+    double ag = atan2(sina, cosa);
     return ag;
   }
-  double angle_deg(const  BDS_Vector &v) const
+  double angle_deg(const BDS_Vector &v) const
   {
     return angle(v) * 180. / M_PI;
   }
-  double operator * (const  BDS_Vector &v) const
+  double operator * (const BDS_Vector &v) const
   {
-    return (x*v.x+y*v.y+z*v.z);
+    return (x * v.x + y * v.y + z * v.z);
   }
-  BDS_Vector(const BDS_Point &p2,const BDS_Point &p1);
+  BDS_Vector(const BDS_Point &p2, const BDS_Point &p1);
   
   BDS_Vector(const double X=0., const double Y=0., const double Z=0.)
-    : x(X),y(Y),z(Z)
+    : x(X), y(Y), z(Z)
   {
   }
   static double t;
@@ -173,16 +173,16 @@ class BDS_Point
   // and is propagated
   double _lcBGM, _lcPTS;
 public:
-  double X,Y,Z; // Real COORDINATES
-  double u,v;   // Parametric COORDINATES
+  double X, Y, Z;
+  double u, v;
   bool config_modified;
   int iD;
   BDS_GeomEntity *g;
   std::list<BDS_Edge*> edges;
 
   // just a transition
-  double & lcBGM  () {return _lcBGM;}
-  double & lc     () {return _lcPTS;}
+  double &lcBGM() { return _lcBGM; }
+  double &lc() { return _lcPTS; }
   
   inline bool operator < (const BDS_Point & other) const
   {
@@ -190,9 +190,9 @@ public:
   }
   inline void del(BDS_Edge *e)
   {
-    std::list<BDS_Edge*>::iterator it  = edges.begin();
+    std::list<BDS_Edge*>::iterator it = edges.begin();
     std::list<BDS_Edge*>::iterator ite = edges.end();
-    while(it!=ite){
+    while(it != ite){
       if(*it == e){
 	edges.erase(it);
 	break;
@@ -200,9 +200,10 @@ public:
       ++it;
     }
   }
-  void getTriangles(std::list<BDS_Face *> &t) const; 	
+  void getTriangles(std::list<BDS_Face *> &t) const;
   BDS_Point(int id, double x=0, double y=0, double z=0)
-    : _lcBGM(1.e22),_lcPTS(1.e22),X(x),Y(y),Z(z),u(0),v(0),config_modified(true),iD(id),g(0)
+    : _lcBGM(1.e22), _lcPTS(1.e22), X(x), Y(y), Z(z), u(0), v(0),
+      config_modified(true), iD(id), g(0)
   {	    
   }
 };
@@ -210,12 +211,12 @@ public:
 class BDS_Edge
 {
   double _length;
-  std::vector <BDS_Face *> _faces;
+  std::vector<BDS_Face*> _faces;
 public:
   bool deleted;
-  BDS_Point *p1,*p2;
+  BDS_Point *p1, *p2;
   BDS_GeomEntity *g;
-  inline BDS_Face* faces(int i) const
+  inline BDS_Face *faces(int i) const
   {
     return _faces [i];
   }
@@ -228,13 +229,13 @@ public:
     return _faces.size();
   }
   int numTriangles() const ;
-  inline BDS_Point * commonvertex(const BDS_Edge *other) const
+  inline BDS_Point *commonvertex(const BDS_Edge *other) const
   {
     if(p1 == other->p1 || p1 == other->p2) return p1;
     if(p2 == other->p1 || p2 == other->p2) return p2;
     return 0;
   }
-  inline BDS_Point * othervertex(const BDS_Point *p) const
+  inline BDS_Point *othervertex(const BDS_Point *p) const
   {
     if(p1 == p) return p2;
     if(p2 == p) return p1;
@@ -244,7 +245,7 @@ public:
   {
     _faces.push_back(f);
   }
-  inline bool operator < (const BDS_Edge & other) const
+  inline bool operator < (const BDS_Edge &other) const
   {
     if(*other.p1 < *p1) return true;
     if(*p1 < *other.p1) return false;
@@ -260,27 +261,30 @@ public:
   }
   inline void del(BDS_Face *t)
   {
-    _faces.erase(std::remove_if(_faces.begin(),_faces.end() , std::bind2nd(std::equal_to<BDS_Face*>(), t)) , 
+    _faces.erase(std::remove_if(_faces.begin(),_faces.end(), 
+				std::bind2nd(std::equal_to<BDS_Face*>(), t)), 
 		 _faces.end());
   }
   
   void oppositeof(BDS_Point * oface[2]) const; 
   
-  void update ()
+  void update()
   {
-    _length = sqrt((p1->X-p2->X)*(p1->X-p2->X)+(p1->Y-p2->Y)*(p1->Y-p2->Y)+(p1->Z-p2->Z)*(p1->Z-p2->Z));
+    _length = sqrt((p1->X - p2->X) * (p1->X - p2->X) + 
+		   (p1->Y - p2->Y) * (p1->Y - p2->Y) + 
+		   (p1->Z - p2->Z) * (p1->Z - p2->Z));
   }
 
   BDS_Edge(BDS_Point *A, BDS_Point *B)
-    : deleted(false),g(0)
+    : deleted(false), g(0)
   {	    
     if(*A < *B){
-      p1=A;
-      p2=B;
+      p1 = A;
+      p2 = B;
     }
     else{
-      p1=B;
-      p2=A;
+      p1 = B;
+      p2 = A;
     }
     p1->edges.push_back(this);
     p2->edges.push_back(this);
@@ -288,47 +292,43 @@ public:
   }
 };
 
-
 class BDS_Face
 {
 public:
   bool deleted;
-  BDS_Edge *e1,*e2,*e3,*e4;
+  BDS_Edge *e1, *e2, *e3, *e4;
   BDS_GeomEntity *g;
-  inline int numEdges () const {return e4?4:3;}
+  inline int numEdges () const { return e4 ? 4 : 3; }
   inline void getNodes(BDS_Point *n[4]) const
   {
-    if (!e4)
-      {
-	n[0] = e1->commonvertex(e3);
-	n[1] = e1->commonvertex(e2);
-	n[2] = e2->commonvertex(e3);
-	n[3] = 0;
-      }
-    else
-      {
-	n[0] = e1->commonvertex(e4);
-	n[1] = e1->commonvertex(e2);
-	n[2] = e2->commonvertex(e3);
-	n[3] = e3->commonvertex(e4);
-      }
+    if (!e4){
+      n[0] = e1->commonvertex(e3);
+      n[1] = e1->commonvertex(e2);
+      n[2] = e2->commonvertex(e3);
+      n[3] = 0;
+    }
+    else{
+      n[0] = e1->commonvertex(e4);
+      n[1] = e1->commonvertex(e2);
+      n[2] = e2->commonvertex(e3);
+      n[3] = e3->commonvertex(e4);
+    }
   }
   
   BDS_Face(BDS_Edge *A, BDS_Edge *B, BDS_Edge *C,BDS_Edge *D = 0)
-    : deleted(false) ,e1(A),e2(B),e3(C),e4(D),g(0)
+    : deleted(false), e1(A), e2(B), e3(C), e4(D), g(0)
   {	
     e1->addface(this);
     e2->addface(this);
     e3->addface(this);
-    if(e4)e4->addface(this);
+    if(e4) e4->addface(this);
   }
 };
 
-
 class GeomLessThan
 {
 public:
-  bool operator()(const BDS_GeomEntity* ent1, const BDS_GeomEntity* ent2) const
+  bool operator()(const BDS_GeomEntity *ent1, const BDS_GeomEntity *ent2) const
   {
     return *ent1 < *ent2;
   }
@@ -337,7 +337,7 @@ public:
 class PointLessThan
 {
 public:
-  bool operator()(const BDS_Point* ent1, const BDS_Point* ent2) const
+  bool operator()(const BDS_Point *ent1, const BDS_Point *ent2) const
   {
     return *ent1 < *ent2;
   }
@@ -347,7 +347,7 @@ class PointLessThanLexicographic
 {
 public:
   static double t;
-  bool operator()(const BDS_Point* ent1, const BDS_Point* ent2) const
+  bool operator()(const BDS_Point *ent1, const BDS_Point *ent2) const
   {
     if(ent1->X - ent2->X  >  t) return true;
     if(ent1->X - ent2->X  < -t) return false;
@@ -361,22 +361,21 @@ public:
 class EdgeLessThan
 {
 public:
-  bool operator()(const BDS_Edge* ent1, const BDS_Edge* ent2) const
+  bool operator()(const BDS_Edge *ent1, const BDS_Edge *ent2) const
   {
     return *ent1 < *ent2;
   }
 };
 
-
 class BDS_SwapEdgeTest
 {
  public:
-  virtual bool operator() (BDS_Point *p1,BDS_Point *p2,
-			   BDS_Point *q1,BDS_Point *q2) const = 0; 
-  virtual bool operator() (BDS_Point *p1,BDS_Point *p2, BDS_Point *p3,
-			   BDS_Point *q1,BDS_Point *q2, BDS_Point *q3,
-			   BDS_Point *op1,BDS_Point *op2, BDS_Point *op3,
-			   BDS_Point *oq1,BDS_Point *oq2, BDS_Point *oq3) const = 0; 
+  virtual bool operator() (BDS_Point *p1, BDS_Point *p2,
+			   BDS_Point *q1, BDS_Point *q2) const = 0; 
+  virtual bool operator() (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
+			   BDS_Point *q1, BDS_Point *q2, BDS_Point *q3,
+			   BDS_Point *op1, BDS_Point *op2, BDS_Point *op3,
+			   BDS_Point *oq1, BDS_Point *oq2, BDS_Point *oq3) const = 0; 
   virtual ~BDS_SwapEdgeTest(){}
 };
 
@@ -384,13 +383,13 @@ class BDS_SwapEdgeTestQuality : public BDS_SwapEdgeTest
 {
   bool testQuality, testSmallTriangles;
  public:
-  BDS_SwapEdgeTestQuality (bool a, bool b=true) : testQuality(a),testSmallTriangles(b){}
-  virtual bool operator() (BDS_Point *p1,BDS_Point *p2,
-			   BDS_Point *q1,BDS_Point *q2) const ; 
-  virtual bool operator() (BDS_Point *p1,BDS_Point *p2, BDS_Point *p3,
-			   BDS_Point *q1,BDS_Point *q2, BDS_Point *q3,
-			   BDS_Point *op1,BDS_Point *op2, BDS_Point *op3,
-			   BDS_Point *oq1,BDS_Point *oq2, BDS_Point *oq3) const ; 
+  BDS_SwapEdgeTestQuality(bool a, bool b=true) : testQuality(a), testSmallTriangles(b) {}
+  virtual bool operator() (BDS_Point *p1, BDS_Point *p2,
+			   BDS_Point *q1, BDS_Point *q2) const ; 
+  virtual bool operator() (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
+			   BDS_Point *q1, BDS_Point *q2, BDS_Point *q3,
+			   BDS_Point *op1, BDS_Point *op2, BDS_Point *op3,
+			   BDS_Point *oq1, BDS_Point *oq2, BDS_Point *oq3) const ; 
   virtual ~BDS_SwapEdgeTestQuality(){}
 };
 
@@ -398,24 +397,22 @@ struct EdgeToRecover
 {
   int p1,p2;
   GEdge *ge;
-  EdgeToRecover ( int  _p1 ,  int  _p2 , GEdge *_ge) : ge(_ge)
+  EdgeToRecover(int _p1, int _p2, GEdge *_ge) : ge(_ge)
   {
-    if (_p1 < _p2 )
-      {
-	p1 = _p1 ;
-	p2 = _p2 ;
-      }
-    else
-      {
-	p2 = _p1 ;
-	p1 = _p2 ;
-      }
+    if(_p1 < _p2 ){
+      p1 = _p1;
+      p2 = _p2;
+    }
+    else{
+      p2 = _p1;
+      p1 = _p2;
+    }
   }
-  bool operator < ( const EdgeToRecover & other) const
+  bool operator < (const EdgeToRecover &other) const
   {    
-    if ( p1 < other.p1 ) return true;
-    if ( p1 > other.p1 ) return false;
-    if ( p2 < other.p2 ) return true;
+    if(p1 < other.p1) return true;
+    if(p1 > other.p1) return false;
+    if(p2 < other.p2) return true;
     return false;
   }
 };
@@ -423,61 +420,62 @@ struct EdgeToRecover
 class BDS_Mesh 
 {    
 public:
-  int MAXPOINTNUMBER,SNAP_SUCCESS,SNAP_FAILURE;
-  double Min[3],Max[3],LC;
+  int MAXPOINTNUMBER, SNAP_SUCCESS, SNAP_FAILURE;
+  double Min[3], Max[3], LC;
   double scalingU, scalingV;
   BDS_Mesh(int _MAXX = 0) :  MAXPOINTNUMBER(_MAXX),scalingU(1),scalingV(1){}
-  void load(GVertex   *gv); // load in BDS all the meshes of the vertex 
-  void load(GEdge     *ge); // load in BDS all the meshes of the edge 
-  void load(GFace     *gf); // load in BDS all the meshes of the surface 
+  void load(GVertex *gv); // load in BDS all the meshes of the vertex 
+  void load(GEdge *ge); // load in BDS all the meshes of the edge 
+  void load(GFace *gf); // load in BDS all the meshes of the surface 
   virtual ~BDS_Mesh();
   BDS_Mesh(const BDS_Mesh &other);
-  std::set<BDS_GeomEntity*,GeomLessThan> geom; 
-  std::set<BDS_Point*,PointLessThan>     points; 
-  std::list<BDS_Edge*>      edges; 
-  std::list<BDS_Face*>   triangles; 
+  std::set<BDS_GeomEntity*, GeomLessThan> geom; 
+  std::set<BDS_Point*, PointLessThan> points; 
+  std::list<BDS_Edge*> edges; 
+  std::list<BDS_Face*> triangles; 
   // Points
-  BDS_Point * add_point(int num , double x, double y,double z);
-  BDS_Point * add_point(int num , double u, double v , GFace *gf);
+  BDS_Point *add_point(int num, double x, double y, double z);
+  BDS_Point *add_point(int num, double u, double v, GFace *gf);
   void del_point(BDS_Point *p);
   BDS_Point *find_point(int num);
   // Edges
-  BDS_Edge  * add_edge(int p1, int p2);
+  BDS_Edge *add_edge(int p1, int p2);
   void del_edge(BDS_Edge *e);
-  BDS_Edge  *find_edge(int p1, int p2);
-  BDS_Edge  *find_edge(BDS_Point*p1, BDS_Point *p2);
-  BDS_Edge  *find_edge(BDS_Point*p1, int p2);
-  BDS_Edge  *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t)const;
+  BDS_Edge *find_edge(int p1, int p2);
+  BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2);
+  BDS_Edge *find_edge(BDS_Point *p1, int p2);
+  BDS_Edge *find_edge(BDS_Point *p1, BDS_Point *p2, BDS_Face *t) const;
   // Triangles & Quadrangles
-  BDS_Face *add_triangle(int p1, int p2, int p3); 
-  BDS_Face *add_quadrangle(int p1, int p2, int p3,int p4); 
-  BDS_Face *add_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
-  BDS_Face *add_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
+  BDS_Face *add_triangle(int p1, int p2, int p3);
+  BDS_Face *add_quadrangle(int p1, int p2, int p3, int p4); 
+  BDS_Face *add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3);
+  BDS_Face *add_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
   void del_face(BDS_Face *t);
-  BDS_Face  *find_triangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3);
-  BDS_Face  *find_quadrangle(BDS_Edge *e1,BDS_Edge *e2,BDS_Edge *e3,BDS_Edge *e4);
+  BDS_Face *find_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3);
+  BDS_Face *find_quadrangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3, BDS_Edge *e4);
   // Geom entities
   void add_geom(int degree, int tag);
   BDS_GeomEntity *get_geom(int p1, int p2);
   // 2D operators
-  BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0, std::set<EdgeToRecover> *not_recovered = 0);
-  bool swap_edge(BDS_Edge *, const BDS_SwapEdgeTest &theTest);
-  bool collapse_edge_parametric(BDS_Edge *, BDS_Point*);
-  void snap_point(BDS_Point* , BDS_Mesh *geom = 0);
-  bool smooth_point(BDS_Point* , BDS_Mesh *geom = 0);
-  bool smooth_point_parametric(BDS_Point * p, GFace *gf);
-  bool smooth_point_centroid(BDS_Point * p, GFace *gf, bool test_quality = false);
+  BDS_Edge *recover_edge(int p1, int p2, std::set<EdgeToRecover> *e2r=0,
+			 std::set<EdgeToRecover> *not_recovered=0);
+  bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest);
+  bool collapse_edge_parametric(BDS_Edge*, BDS_Point*);
+  void snap_point(BDS_Point*, BDS_Mesh *geom = 0);
+  bool smooth_point(BDS_Point*, BDS_Mesh *geom = 0);
+  bool smooth_point_parametric(BDS_Point *p, GFace *gf);
+  bool smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality=false);
   bool move_point(BDS_Point *p , double X, double Y, double Z);
   bool split_edge(BDS_Edge *, BDS_Point *);
   bool split_face(BDS_Face *, BDS_Point *);
-  bool edge_constraint    ( BDS_Point *p1, BDS_Point *p2 );
-  bool recombine_edge    ( BDS_Edge *e );
+  bool edge_constraint(BDS_Point *p1, BDS_Point *p2);
+  bool recombine_edge(BDS_Edge *e);
   // Global operators
   void cleanup();
-  void recombineIntoQuads (const double angle, GFace *gf);
+  void recombineIntoQuads(const double angle, GFace *gf);
 };
 
-void outputScalarField(std::list < BDS_Face * >t, const char *fn, int param, GFace *gf = 0);
-void recur_tag(BDS_Face * t, BDS_GeomEntity * g);
+void outputScalarField(std::list<BDS_Face*> t, const char *fn, int param, GFace *gf=0);
+void recur_tag(BDS_Face *t, BDS_GeomEntity *g);
 
 #endif