diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index fbdb4705b04aaf66f33d73713cb616614b52fe38..95c98169325c98deff8278833971a00a630ecf13 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -58,7 +58,7 @@ static void class_selectgface_cb(Fl_Widget *w, void *data)
     Draw();
 
     Msg::StatusBar(3, false, "Select Model Face\n"
-        "[Press 'e' to end selection or 'q' to abort]");
+		   "[Press 'e' to end selection or 'q' to abort]");
     
     char ib = GUI::instance()->selectEntity(ENT_SURFACE);
     if(ib == 'l') {
@@ -98,7 +98,7 @@ static void class_deleteedge_cb(Fl_Widget *w, void *data)
     Draw();
 
     Msg::StatusBar(3, false, "Select Elements\n"
-        "[Press 'e' to end selection or 'q' to abort]");
+		   "[Press 'e' to end selection or 'q' to abort]");
     
     char ib = GUI::instance()->selectEntity(ENT_ALL);
     if(ib == 'l') {
@@ -275,6 +275,7 @@ static GEdge *getNewModelEdge(GFace *gf1, GFace *gf2,
     newEdges.find(std::make_pair<int, int>(i1, i2));
   if(it == newEdges.end()){
     discreteEdge *temporary = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+    printf("add new edge gf1=%d gf2=%d \n", t1, t2);
     GModel::current()->add(temporary);
     newEdges[std::make_pair<int, int>(i1, i2)] = temporary;
     return temporary;
@@ -290,6 +291,7 @@ static void recurClassifyEdges(MTri3 *t,
                                std::map<std::pair<int, int>, GEdge*> &newEdges)
 {
   if(!t->isDeleted()){
+
     t->setDeleted(true);
     GFace *gf1 = reverse[t->tri()];
     for(int i = 0; i < 3; i++){
@@ -310,6 +312,7 @@ static void recurClassifyEdges(MTri3 *t,
       if(tn)
         recurClassifyEdges(tn, reverse, lines, touched, newEdges);
     }
+
   }
 }
 
@@ -367,11 +370,97 @@ static void class_color_cb(Fl_Widget* w, void* data)
     
     it = tris.begin();
     
+    //classify edges that are bound by different GFaces
+    //--------------------------------------------------
     std::map<std::pair<int, int>, GEdge*> newEdges;
     std::set<MLine*> touched;
     recurClassifyEdges(*it, reverse, lines, touched, newEdges);
     GModel::current()->remove(e->saved);
 
+    //check if new edges should not be splitted 
+    //splitted if composed of several open or closed edges
+    //-----------------------------------------------------
+    for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it){
+
+      GEdge *ge = it->second;
+      printf("new edge with tag %d \n", ge->tag());
+
+      std::list<MLine*> segments;
+      for (int i=0; i < ge->lines.size(); i++){
+	segments.push_back(ge->lines[i]);
+      }
+
+      //for each actual GEdge
+      while (! segments.empty()) {
+	std::vector<MLine*> myLines;
+	std::list<MLine*>::iterator it = segments.begin();
+
+	MVertex *vB = (*it)->getVertex(0);
+	MVertex *vE = (*it)->getVertex(1);
+	myLines.push_back(*it);
+	segments.erase(it);
+	it++;
+
+	//printf("***candidate mline %d %d of size \n", vB->getNum(), vE->getNum(), segments.size());
+
+   	for (int i=0; i<2; i++) {
+
+	  bool found = false;
+	  for (std::list<MLine*>::iterator it = segments.begin() ; it != segments.end(); ++it){	
+	    MVertex *v1 = (*it)->getVertex(0);
+	    MVertex *v2 = (*it)->getVertex(1);
+	    //printf("mline %d %d \n", v1->getNum(), v2->getNum());
+
+	    if ( v1 == vE  ){
+	      //printf("->push back this mline \n");
+	      myLines.push_back(*it);
+	      segments.erase(it);
+	      vE = v2;
+	      i = -1;
+	    }
+	    else if ( v2 == vE){
+	      //printf("->push back this mline \n");
+	      myLines.push_back(*it);
+	      segments.erase(it);
+	      vE = v1;
+	      i=-1;
+	    }
+	  }
+
+	  if (vB == vE) break;
+
+	  if (segments.empty()) break;
+
+	  //printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
+	  MVertex *temp = vB;
+	  vB = vE;
+	  vE = temp;
+	  //printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
+
+	}
+	
+// 	printf("************ CANDIDATE NEW EDGE \n");
+// 	for (std::vector<MLine*>::iterator it = myLines.begin() ; it != myLines.end() ; ++it){
+// 	  MVertex *v1 = (*it)->getVertex(0);
+// 	  MVertex *v2 = (*it)->getVertex(1);
+// 	  printf("Line %d %d \n", v1->getNum(), v2->getNum());
+// 	}
+	GEdge *newGe = new discreteEdge(GModel::current(), maxEdgeNum() + 1, 0, 0);
+	newGe->lines.insert(newGe->lines.end(), myLines.begin(), myLines.end());
+	GModel::current()->add(newGe);
+	//printf("create new edge with tag =%d\n", maxEdgeNum());
+	
+      }//end for each actual GEdge
+  
+    }
+
+    for (std::map<std::pair<int, int>, GEdge*>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it){
+      GEdge *ge = it->second;
+      GModel::current()->remove(ge);
+    }
+
+
+
     while(it != tris.end()){
       delete *it;
       ++it;
@@ -451,7 +540,7 @@ static void class_select_cb(Fl_Widget *w, void *data)
     Draw();
 
     Msg::StatusBar(3, false, "Select Elements\n"
-        "[Press 'e' to end selection or 'q' to abort]");
+		   "[Press 'e' to end selection or 'q' to abort]");
     
     char ib = GUI::instance()->selectEntity(ENT_ALL);
     if(ib == 'l') {
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index d47cf368577e42be37bb9e2e37bfdade73d3dca1..3eb3687f4e0f936878a3d24550b0ae03a5f165a0 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -92,22 +92,26 @@ void GFaceCompound::parametrize() const
 
 void GFaceCompound::getBoundingEdges()
 {
-  //printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
+  printf("***** In GFaceCompound: size U0=%d, v0=%d\n ", _U0.size(), _V0.size());
 
+  //in case the bounding edges are explicitely given
   if (_U0.size()){
     std::list<GEdge*> :: const_iterator it = _U0.begin();
     for ( ; it != _U0.end() ; ++it){
       l_edges.push_back(*it);
+      //printf("U0 for edge %d, add face %d \n", (*it)->tag(), this->tag());
       (*it)->addFace(this);
     }
     it = _V0.begin();
     for ( ; it != _V0.end() ; ++it){
       l_edges.push_back(*it);
+      //printf("V0 for edge %d, add face %d \n", (*it)->tag(), this->tag());
       (*it)->addFace(this);
     }
     return;
   }
 
+  // in case the bounding edges are not given Boundary { {} };
   std::set<GEdge*> _unique;
   std::multiset<GEdge*> _touched;
   std::list<GFace*>::iterator it = _compound.begin();
@@ -124,8 +128,7 @@ void GFaceCompound::getBoundingEdges()
     std::list<GEdge*> :: iterator ite = ed.begin();
     for ( ; ite != ed.end() ; ++ite){
       if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1) {	
-	_unique.insert(*ite);
-      }
+	_unique.insert(*ite);      }
     }    
   }    
 
@@ -136,7 +139,34 @@ void GFaceCompound::getBoundingEdges()
     (*itf)->addFace(this);
   }
 
-  _U0 = l_edges;
+  //find a closed loop for assigning boundary conditions
+  std::list<GEdge*> _loop;
+  std::map<GVertex*,GEdge*> tempv;
+  std::set<GEdge*>::iterator its = _unique.begin();
+  GVertex *vB = (*its)->getBeginVertex();
+  GVertex *v = (*its)->getEndVertex();
+  _loop.push_back(*its);
+  printf("boundary add edge=%d \n", (*its)->tag());
+  //printf("edge=%d  vB=%d v =%d\n", (*its)->tag(), vB->tag(), v->tag());
+  its++;
+  for ( ; its != _unique.end() ; ++its){
+    GVertex *v1 = (*its)->getBeginVertex();
+    GVertex *v2 = (*its)->getEndVertex();
+    if (v1 == v) {
+      printf("boundary add edge=%d \n", (*its)->tag());
+      _loop.push_back(*its);
+      v = v2;
+    }
+    else if (v2 == v) {
+      printf("boundaryn add edge=%d \n", (*its)->tag());
+      _loop.push_back(*its);
+      v = v1;
+    }
+    if (v == vB) break;
+  }
+ 
+  _U0 = _loop;
+
 }
 
 GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
@@ -148,6 +178,11 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
   if (!_U0.size()) _type = UNITCIRCLE;
   else if (!_V1.size()) _type = UNITCIRCLE;
   else _type = SQUARE;
+
+  for (std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+    if (!(*it)) Msg::Error("Incorrect face in compound surface %d\n", tag);
+  }
+
 }
 
 GFaceCompound::~GFaceCompound()
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index dbd02253abdf3ba8530bae7549ce778517ecdadc..b638d242ea5d0ea0f3f2c2c6f3ef7626d249cb63 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -113,6 +113,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
   if(part) m->getMeshPartitions().insert(part);
 }
 
+
 int GModel::readMSH(const std::string &name)
 {
   FILE *fp = fopen(name.c_str(), "rb");
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index bf40406911937f28a172d0ffad54c9ab8b81b40f..21b5212dfba36f6a35d81e65c491b0d026994ddb 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -146,10 +146,9 @@ void discreteEdge::setBoundVertices()
     for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); iter != boundv.end(); iter++){
       MVertex* vE = (iter)->first;
       bool existVertex  = false;
-      //for (std::vector<discreteVertex*>::iterator point = vertices.begin(); point != vertices.end(); point++) {
       for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
-	//printf("Discrete point =%d bound vE=%d\n", (*point)->tag(), vE->getNum());
-	if ((*point)->tag() == vE->getNum()){
+	//printf("Discrete point =%d %d  bound vE=%d %d\n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum(), vE->getIndex());
+	if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
 	  bound_vertices.push_back((*point));
 	  existVertex = true;
 	  break;
@@ -174,7 +173,7 @@ void discreteEdge::setBoundVertices()
     v1 = bound_vertices[1];
   }
   else if (boundv.size()==0){ 
-    //printf("Found a closed Curve \n");
+    //printf("Found a closed Curve add vertex \n");
 
     std::vector<MLine*>::const_iterator it = lines.begin();
     MVertex* vE = (*it)->getVertex(0);
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 439a7ec1bb18b2df8bb3c59b556db882fa8fd051..9f2c213eb831ad3801897754d299e28867fe1973 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -22,9 +22,11 @@ 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
+
   std::list<MEdge> bound_edges;
   for (int iFace = 0; iFace  < getNumMeshElements() ; iFace++) {
     std::vector<MVertex*> fv;
@@ -63,16 +65,6 @@ void discreteFace::setBoundEdges(std::vector<int> tagEdges)
 
   // printf("***** In discrete Face:  \n");
 
-  
-//   for(GModel::eiter it = model()->firstEdge(); it != model()->lastEdge(); ++it){
-// //   for (std::vector<discreteEdge*>::iterator it = bound_edges.begin(); it != bound_edges.end(); it++) {
-//     l_edges.push_back(*it);
-//     l_dirs.push_back(1);
-//     (*it)->addFace(this);
-//     printf("add Face %d for edge %d\n",this->tag(), (*it)->tag() );
-//   }
-
-
  for (std::vector<int>::iterator it = tagEdges.begin(); it != tagEdges.end(); it++) {
    GEdge *ge = GModel::current()->getEdgeByTag(abs(*it));
    l_edges.push_back(ge);
diff --git a/benchmarks/2d/xytouv.geo b/benchmarks/2d/xytouv.geo
index 1a5c61cf4bd3aae9e86cc3808e63d7cbf43c8fa4..73873eea5d941a6a8b0eaeef5b22e5dd4f8a3e23 100644
--- a/benchmarks/2d/xytouv.geo
+++ b/benchmarks/2d/xytouv.geo
@@ -1,4 +1,4 @@
-lc = 0.2;
+lc = 0.3;
 Point(1)={0,0,0,lc};
 Point(2)={1,0,0,lc};
 Point(3)={0,1,0,lc};