diff --git a/Geo/ExtractContour.cpp b/Geo/ExtractContour.cpp
index 9ce767b8802ff14b6b03eaf8330965394c3cc65d..ac1f22ec10719fff61c6b19e18497f23e165562f 100644
--- a/Geo/ExtractContour.cpp
+++ b/Geo/ExtractContour.cpp
@@ -1,4 +1,4 @@
-// $Id: ExtractContour.cpp,v 1.3 2004-05-19 04:03:51 geuzaine Exp $
+// $Id: ExtractContour.cpp,v 1.4 2004-06-30 00:57:50 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -24,6 +24,7 @@
 #include "GeoUtils.h"
 #include "CAD.h"
 #include "Mesh.h"
+#include "Numeric.h"
 
 // Note: we use List_ISearchSeq so that the input lists don't get
 // sorted: it's less efficient, but it allows us to do multi-level
@@ -114,6 +115,9 @@ void createEdgeLinks(Tree_T *links)
 
 void orientAndSortEdges(List_T *edges, Tree_T *links)
 {
+  // this orients all the edges in a line loop in a consistent manner
+  // (left- or right-oriented, depending on the orientation of the
+  // first edge) *and* sorts them so that they form a path
   int num;
   lnk lk;
   nxa na;
@@ -185,11 +189,12 @@ int allEdgesLinked(int ed, List_T * edges)
 
   int found = 0;
 
-  if(!Tree_Nbr(points))
+  if(!Tree_Nbr(points)){
     found = 1;
-
-  if(found)
+    // we can only orient things now since we allow to select
+    // disconnected parts of the loop interactively
     orientAndSortEdges(edges, links);
+  }
 
   Tree_Delete(links);
   Tree_Delete(points);
@@ -258,6 +263,60 @@ void createFaceLinks(Tree_T * links)
   List_Delete(temp);
 }
 
+void recurOrientFace(int face, List_T *faces, List_T *available, Tree_T *links)
+{
+  Surface *s = FindSurface(abs(face), THEM);
+  int ori = sign(face);
+
+  for(int i = 0; i < List_Nbr(s->Generatrices); i++){
+    Curve *c;
+    List_Read(s->Generatrices, i, &c);
+    lnk lk;
+    lk.n = abs(c->Num);
+    Tree_Query(links, &lk);
+    for(int j = 0; j < List_Nbr(lk.l); j++){
+      nxa na;
+      List_Read(lk.l, j, &na);
+      int num = abs(na.a);
+      if(num != abs(s->Num) && List_Search(available, &num, fcmp_absint) &&
+	 List_ISearchSeq(faces, &num, fcmp_absint) < 0){
+	Surface *s2 = FindSurface(num, THEM);
+	for(int k = 0; k < List_Nbr(s2->Generatrices); k++){
+	  Curve *c2;
+	  List_Read(s2->Generatrices, k, &c2);
+	  if(abs(c->Num) == abs(c2->Num)){
+	    if(c->Num * c2->Num > 0)
+	      num *= -ori;
+	    else
+	      num *= ori;
+	    List_Add(faces, &num);
+	    recurOrientFace(num, faces, available, links);
+	    break;
+	  }
+	}
+      }
+    }
+  }
+}
+
+void orientFaces(List_T *faces, Tree_T *links)
+{
+  // this orients all the faces in a surface loop in a consistent
+  // manner (all normals pointing inside or outside--depending on the
+  // orientation of the first face)
+
+  List_T *temp = List_Create(List_Nbr(faces), 1, sizeof(int));
+  List_Copy(faces, temp);
+  List_Reset(faces);
+
+  int num;
+  List_Read(temp, 0, &num);
+  List_Add(faces, &num);
+  recurOrientFace(num, faces, temp, links);
+
+  List_Delete(temp);
+}
+
 int allFacesLinked(int fac, List_T * faces)
 {
   Tree_T *links = Tree_Create(sizeof(lnk), complink);
@@ -284,13 +343,21 @@ int allFacesLinked(int fac, List_T * faces)
 
   if(List_ISearchSeq(faces, &fac, fcmp_absint) < 0){
     List_Add(faces, &fac);
+    // Warning: this is correct ONLY if the surfaces are defined
+    // correctly, i.e., if the surface hole boundaries are oriented
+    // consistently with the surface exterior boundaries. There is
+    // currently *nothing* in the code that checks this.
     recurFindLinkedFaces(fac, faces, edges, links);
   }
 
   int found = 0;
 
-  if(!Tree_Nbr(edges))
+  if(!Tree_Nbr(edges)){
     found = 1;
+    // we can only orient things now since we allow to select
+    // disconnected parts of the loop interactively
+    orientFaces(faces, links);
+  }
 
   Tree_Delete(links);
   Tree_Delete(edges);
diff --git a/Geo/GeoUtils.cpp b/Geo/GeoUtils.cpp
index a7fc3e7095767cc608d7e1c2c27e90e759348ada..3645d8809313518d2a68987f059a85aa5cd06b5d 100644
--- a/Geo/GeoUtils.cpp
+++ b/Geo/GeoUtils.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoUtils.cpp,v 1.3 2004-06-26 17:58:14 geuzaine Exp $
+// $Id: GeoUtils.cpp,v 1.4 2004-06-30 00:57:50 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -23,12 +23,14 @@
 #include "Geo.h"
 #include "CAD.h"
 #include "Mesh.h"
+#include "Numeric.h"
 
 extern Mesh *THEM;
 
-// A small function to sort the edges in an EdgeLoop. Without this
-// sort, it is very difficult to write general scriptable surface
-// generation in complex cases
+// This function sorts the edges in an EdgeLoop and detects any
+// subloops. Warning: the input edges are supposed to be *oriented*
+// (Without this sort, it is very difficult to write general
+// scriptable surface generation in complex cases)
 
 void sortEdgesInLoop(int num, List_T *edges)
 {
@@ -94,17 +96,33 @@ void setSurfaceGeneratrices(Surface *s, List_T *loops)
       return;
     }
     else {
-      for(int j = 0; j < List_Nbr(el->Curves); j++) {
-	int ic;
-        List_Read(el->Curves, j, &ic);
-	Curve *c;
-        if(!(c = FindCurve(ic, THEM))) {
-          Msg(GERROR, "Unknown Curve %d", ic);
-          List_Delete(s->Generatrices);
-          return;
-        }
-        else
-          List_Add(s->Generatrices, &c);
+      int ic;
+      Curve *c;
+      if(i == 0){ // exterior boundary
+	for(int j = 0; j < List_Nbr(el->Curves); j++) {
+	  List_Read(el->Curves, j, &ic);
+	  if(!(c = FindCurve(ic, THEM))) {
+	    Msg(GERROR, "Unknown Curve %d", ic);
+	    List_Delete(s->Generatrices);
+	    return;
+	  }
+	  else
+	    List_Add(s->Generatrices, &c);
+	}
+      }
+      else{ // holes, assuming that their orientation is consistent
+	    // with the orientation of the exterior boundary!
+	for(int j = List_Nbr(el->Curves)-1; j >= 0; j--) {
+	  List_Read(el->Curves, j, &ic);
+	  ic *= -1;
+	  if(!(c = FindCurve(ic, THEM))) {
+	    Msg(GERROR, "Unknown Curve %d", ic);
+	    List_Delete(s->Generatrices);
+	    return;
+	  }
+	  else
+	    List_Add(s->Generatrices, &c);
+	}
       }
     }
   }
@@ -140,7 +158,8 @@ void setVolumeSurfaces(Volume *v, List_T * loops)
         }
         else{
           List_Add(v->Surfaces, &s);
-	  int tmp = (is > 0) ? 1 : -1;
+	  int tmp = sign(is);
+	  if(i > 0) tmp *= -1; // this is a hole
 	  List_Add(v->SurfacesOrientations, &tmp);
 	}
       }
diff --git a/Mesh/2D_Mesh.cpp b/Mesh/2D_Mesh.cpp
index 0000cf8d88d27534a8051f13c65ef53853ec60d2..69862d28a4ae43ab786f6717fdb8727a6493dacd 100644
--- a/Mesh/2D_Mesh.cpp
+++ b/Mesh/2D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh.cpp,v 1.66 2004-06-23 22:25:01 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.67 2004-06-30 00:57:50 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -805,6 +805,9 @@ void ReOrientSurfaceMesh(Surface *s)
   if(List_Nbr(c->Vertices) < 2)
     return;
 
+  // take the first edge (even if there are holes, we know that this
+  // edge belongs to the exterior loop, due to the ordering of
+  // s->Generatrices)
   Vertex *v1, *v2;
   List_Read(c->Vertices, 0, &v1);
   List_Read(c->Vertices, 1, &v2);
diff --git a/TODO b/TODO
index 8f3392332eb28476e2c8ae1388f5ff9e610e0f37..41dcc292e4a859d15b63ade3e46f2c0e8bb25cd9 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,9 @@
-$Id: TODO,v 1.52 2004-06-24 07:13:18 geuzaine Exp $
+$Id: TODO,v 1.53 2004-06-30 00:57:49 geuzaine Exp $
+
+add an interactive way to choose the orientation of surfaces in
+surface loops and lines in line loops?
+
+********************************************************************
 
 add ternary operator and <,>,<=,>=,== tests in MathEval