Replaced my orientation hack with a rigorous algorithm... Seems to
work nicely (and is general).

Graphics/Geom.cpp

-// $Id: Geom.cpp,v 1.67 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: Geom.cpp,v 1.68 2004-06-23 18:52:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -207,7 +207,7 @@ void Draw_Curve(void *a, void *b)
 
 // Surfaces
 
-void put_Z(Vertex * v, Surface * s)
+void putZ(Vertex * v, Surface * s)
 {
   Vertex V;
   V.Pos.X = s->a;
@@ -255,6 +255,45 @@ int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z,
   return 0;
 }
 
+void getPlaneSurfaceNormal(Surface *s, double x, double y, double z, double n[3])
+{
+  double t1[3], t2[3];
+  Curve *c;
+
+  for(int i = 0; i < 3; i++)
+    n[i] = 0.0;
+
+  if(s->Typ != MSH_SURF_PLAN)
+    return;
+
+  // We cannot use s->plan or InterpolateSurface here (they both rely
+  // on the mean plane, which borks the orientation). So we need to
+  // use this trick:
+
+  if(List_Nbr(s->Generatrices) > 1){
+    List_Read(s->Generatrices, 0, &c);
+    t1[0] = x - c->beg->Pos.X;
+    t1[1] = y - c->beg->Pos.Y;
+    t1[2] = z - c->beg->Pos.Z;
+    for(int i = 1; i < List_Nbr(s->Generatrices); i++){
+      List_Read(s->Generatrices, i, &c);
+      t2[0] = x - c->beg->Pos.X;
+      t2[1] = y - c->beg->Pos.Y;
+      t2[2] = z - c->beg->Pos.Z;
+      prodve(t1, t2, n);
+      if(norme(n))
+	break;
+    }
+  }
+  if(List_Nbr(s->Generatrices) < 2 || !norme(n)){
+    Msg(WARNING, "Reverting to mean plane normal for surface %d", s->Num);
+    n[0] = s->a;
+    n[1] = s->b;
+    n[2] = s->c;
+    norme(n);
+  }
+}
+
 void Draw_Triangulated_Surface(Surface * s)
 {
   int k = 0;
@@ -344,7 +383,7 @@ void Draw_Plane_Surface(Surface * s)
 
     for(i = 0; i < 4; i++) {
       V[i].Pos.Z = 0.0;
-      put_Z(&V[i], s);
+      putZ(&V[i], s);
     }
 
     n[0] = s->plan[2][0];
@@ -441,21 +480,17 @@ void Draw_Plane_Surface(Surface * s)
     if(CTX.geom.normals) {
       List_Read(s->Orientations, 0, &vv1);
       List_Read(s->Orientations, 1, &vv2);
-      // don't rely on MeanPlane: teh orientation is arbotrary
-      //n[0] = s->plan[2][0];
-      //n[1] = s->plan[2][1];
-      //n[2] = s->plan[2][2];
-      //norme(n);
-      Get_SurfaceNormal(s, n);
+      double x = (vv1.Pos.X + vv2.Pos.X) / 2.;
+      double y = (vv1.Pos.Y + vv2.Pos.Y) / 2.;
+      double z = (vv1.Pos.Z + vv2.Pos.Z) / 2.;
+      getPlaneSurfaceNormal(s, x, y, z, n);
       n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
       n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
       n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
       glColor4ubv((GLubyte *) & CTX.color.geom.normals);
       Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
 		  CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius, 
-		  (vv2.Pos.X + vv1.Pos.X) / 2., (vv2.Pos.Y + vv1.Pos.Y) / 2., 
-		  (vv2.Pos.Z + vv1.Pos.Z) / 2., n[0], n[1], n[2], NULL, 
-		  CTX.geom.light);
+		  x, y, z, n[0], n[1], n[2], NULL, CTX.geom.light);
     }
 
   }
@@ -518,16 +553,25 @@ void Draw_NonPlane_Surface(Surface * s)
   }
 
   if(CTX.geom.normals) {
-    Vertex v = InterpolateSurface(s, 0.5, 0.5, 0, 0);
-    double n[3];
-    Get_SurfaceNormal(s, n);
+    Vertex v1 = InterpolateSurface(s, 0.5, 0.5, 0, 0);
+    Vertex v2 = InterpolateSurface(s, 0.6, 0.5, 0, 0);
+    Vertex v3 = InterpolateSurface(s, 0.5, 0.6, 0, 0);
+    double t1[3], t2[3], n[3];
+    t1[0] = v2.Pos.X - v1.Pos.X;
+    t1[1] = v2.Pos.Y - v1.Pos.Y;
+    t1[2] = v2.Pos.Z - v1.Pos.Z;
+    t2[0] = v3.Pos.X - v1.Pos.X;
+    t2[1] = v3.Pos.Y - v1.Pos.Y;
+    t2[2] = v3.Pos.Z - v1.Pos.Z;
+    prodve(t1, t2, n);
+    norme(n);
     n[0] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[0];
     n[1] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[1];
     n[2] *= CTX.geom.normals * CTX.pixel_equiv_x / CTX.s[2];
     glColor4ubv((GLubyte *) & CTX.color.geom.normals);
     Draw_Vector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
 		CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
-		v.Pos.X, v.Pos.Y, v.Pos.Z, n[0], n[1], n[2], NULL,
+		v1.Pos.X, v1.Pos.Y, v1.Pos.Z, n[0], n[1], n[2], NULL,
 		CTX.geom.light);
   }
 }

Mesh/2D_Mesh.cpp

-// $Id: 2D_Mesh.cpp,v 1.64 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.65 2004-06-23 18:52:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -789,6 +789,8 @@ void ActionInvertQua(void *a, void *b)
   q->V[3] = tmp;
 }
 
+int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z, double n[3]);
+
 int isMiddlePointOnPlanarSurface(Surface *s, Vertex *v1, Vertex *v2)
 {
   double n[3] = {s->a, s->b, s->c};
@@ -856,52 +858,38 @@ void Get_SurfaceNormal(Surface *s, double n[3])
   }
 }
 
-// Horrible (tm) hack to orient the elements correctly. This is
-// *definitely* not the best way to do it, but I don't have time to
-// look into the issue right now.
+// A little routine to re-orient the surface meshe to match the
+// geometrical orientation. It would be better to actually generate
+// the meshes correctly in the first place (!), but this is handy
+// since we use many different algorithms since the mean plane
+// computation doesn't take the original orientation into account.
 
 void ReOrientSurfaceMesh(Surface *s)
 {
-  double t1[3], t2[3], n1[3], n2[3], res;
-  Simplex *simp;
-  Quadrangle *quad;
-
-  Get_SurfaceNormal(s, n1);
-
-  if(Tree_Right(s->Simplexes, &simp)){
-    t1[0] = simp->V[1]->Pos.X - simp->V[0]->Pos.X;
-    t1[1] = simp->V[1]->Pos.Y - simp->V[0]->Pos.Y;
-    t1[2] = simp->V[1]->Pos.Z - simp->V[0]->Pos.Z;
-    t2[0] = simp->V[2]->Pos.X - simp->V[0]->Pos.X;
-    t2[1] = simp->V[2]->Pos.Y - simp->V[0]->Pos.Y;
-    t2[2] = simp->V[2]->Pos.Z - simp->V[0]->Pos.Z;
-    prodve(t1, t2, n2);
-    norme(n2);
-    prosca(n1, n2, &res);
-    if(res < 0.0){
-      Msg(DEBUG, "Inverting triangles in %s surface %d (res = %g)",
-	  (s->Typ == MSH_SURF_PLAN) ? "Plane" : "NonPlane", s->Num, res);
+  Curve *c;
+  List_Read(s->Generatrices, 0, &c);
+  Vertex *v1, *v2;
+  List_Read(c->Vertices, 0, &v1);
+  List_Read(c->Vertices, 1, &v2);
+  Edge edge;
+  edge.V[0] = v1;
+  edge.V[1] = v2;
+
+  if(Tree_Nbr(s->Simplexes)){
+    EdgesContainer edges;
+    edges.AddSimplexTree(s->Simplexes);
+    Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
+    if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num)
       Tree_Action(s->Simplexes, ActionInvertTri);
   }
-  }
-
-  if(Tree_Right(s->Quadrangles, &quad)){
-    t1[0] = quad->V[1]->Pos.X - quad->V[0]->Pos.X;
-    t1[1] = quad->V[1]->Pos.Y - quad->V[0]->Pos.Y;
-    t1[2] = quad->V[1]->Pos.Z - quad->V[0]->Pos.Z;
-    t2[0] = quad->V[2]->Pos.X - quad->V[0]->Pos.X;
-    t2[1] = quad->V[2]->Pos.Y - quad->V[0]->Pos.Y;
-    t2[2] = quad->V[2]->Pos.Z - quad->V[0]->Pos.Z;
-    prodve(t1, t2, n2);
-    norme(n2);
-    prosca(n1, n2, &res);
-    if(res < 0.0){
-      Msg(DEBUG, "Inverting quads in %s surface %d (res = %g)",
-	  (s->Typ == MSH_SURF_PLAN) ? "Plane" : "NonPlane", s->Num, res);
+  if(Tree_Nbr(s->Quadrangles)){
+    EdgesContainer edges;
+    edges.AddQuadrangleTree(s->Quadrangles);
+    Edge *e = (Edge*)Tree_PQuery(edges.AllEdges, &edge);
+    if(e && e->V[0]->Num == v2->Num && e->V[1]->Num == v1->Num)
       Tree_Action(s->Quadrangles, ActionInvertQua);
   }
 }
-}
 
 void Maillage_Surface(void *data, void *dum)
 {
@@ -1006,7 +994,6 @@ void Maillage_Surface(void *data, void *dum)
     End_Surface(s->Support, 0);
     End_Surface(s, 0);
     
-    // FIXME: big hack
     ReOrientSurfaceMesh(s);
   }
 

Mesh/3D_Extrude.cpp

-// $Id: 3D_Extrude.cpp,v 1.81 2004-06-23 03:57:43 geuzaine Exp $
+// $Id: 3D_Extrude.cpp,v 1.82 2004-06-23 18:52:45 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -1046,7 +1046,6 @@ int Extrude_Mesh(Surface * s)
   Tree_Action(s->Simplexes, AddSimVertsInSurf);
   Tree_Action(s->Quadrangles, AddQuadVertsInSurf);
 
-  // FIXME: big hack
   ReOrientSurfaceMesh(s);
 
   return true;

Mesh/Mesh.h

@@ -476,8 +476,6 @@ void Projette_Plan_Moyen(void *a, void *b);
 void Projette_Inverse(void *a, void *b);
 void Freeze_Vertex(void *a, void *b);
 void deFreeze_Vertex(void *a, void *b);
-int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z, double n[3]);
-void Get_SurfaceNormal(Surface *s, double n[3]);
 void ReOrientSurfaceMesh(Surface *s);
 
 double Lc_XYZ(double X, double Y, double Z, Mesh * m);

doc/VERSIONS

-$Id: VERSIONS,v 1.226 2004-06-23 04:02:38 geuzaine Exp $
+$Id: VERSIONS,v 1.227 2004-06-23 18:52:45 geuzaine Exp $
 
 New since 1.53: fixed UNV output; make Layers' region numbering
 consistent between lines/surfaces/volumes; fixed home directory
 problem on Win98; new Plugin(CutParametric); the default project file
 is now created in the home directory if no current directory is
 defined (e.g. when double-clicking on the icon on Windows/MacOS);
-(partial) fix of the geometry/mesh surface orientation discrepancies;
-small bug fixes and cleanups.
+fixed the discrepancy between the orientation of geometrical surfaces
+and the associated surface meshes; small bug fixes and cleanups.
 
 New in 1.53: completed support for second order elements in the mesh
 module (lines, triangles, quadrangles, tetrahedra, hexahedra, prisms