Apply the surface re-orientation hack to quadrangles (and to extruded
meshes!)

Graphics/Geom.cpp

-// $Id: Geom.cpp,v 1.66 2004-06-22 17:34:10 geuzaine Exp $
+// $Id: Geom.cpp,v 1.67 2004-06-23 03:57:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -285,8 +285,6 @@ void Draw_Triangulated_Surface(Surface * s)
   }  
 }
 
-void Get_SurfaceNormal(Surface *s, double n[3]);
-
 void Draw_Plane_Surface(Surface * s)
 {
   int i, j, k;

Mesh/2D_Mesh.cpp

-// $Id: 2D_Mesh.cpp,v 1.63 2004-06-22 22:10:11 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.64 2004-06-23 03:57:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -773,23 +773,21 @@ void ActionEndTheCurve(void *a, void *b)
   End_Curve(c);
 }
 
-void ActionInvertTriQua(void *a, void *b)
+void ActionInvertTri(void *a, void *b)
 {
   Simplex *s = *(Simplex **) a;
-  Vertex *tmp;
-  if(s->V[3]){
-    tmp = s->V[1];
-    s->V[1] = s->V[3];
-    s->V[3] = tmp;
-  }
-  else if(s->V[2]){
-    tmp = s->V[1];
+  Vertex *tmp = s->V[1];
   s->V[1] = s->V[2];
   s->V[2] = tmp;
 }
-}
 
-int isPointOnPlanarSurface(Surface * S, double X, double Y, double Z, double n[3]);
+void ActionInvertQua(void *a, void *b)
+{
+  Quadrangle *q = *(Quadrangle **) a;
+  Vertex *tmp = q->V[1];
+  q->V[1] = q->V[3];
+  q->V[3] = tmp;
+}
 
 int isMiddlePointOnPlanarSurface(Surface *s, Vertex *v1, Vertex *v2)
 {
@@ -858,15 +856,19 @@ 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.
+
 void ReOrientSurfaceMesh(Surface *s)
 {
-  // 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.
-  Simplex *simp;
-  if(Tree_Right(s->Simplexes, &simp)){
   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;
@@ -875,19 +877,28 @@ void ReOrientSurfaceMesh(Surface *s)
     t2[2] = simp->V[2]->Pos.Z - simp->V[0]->Pos.Z;
     prodve(t1, t2, n2);
     norme(n2);
-    /*
-      printf("n1=%g %g %g\n", n1[0], n1[1], n1[2]);
-      printf("n2=%g %g %g (elt: (%g,%g,%g) (%g,%g,%g) (%g,%g,%g))\n", 
-	     n2[0], n2[1], n2[2],
-	     simp->V[0]->Pos.X, simp->V[0]->Pos.Y, simp->V[0]->Pos.Z,
-	     simp->V[1]->Pos.X, simp->V[1]->Pos.Y, simp->V[1]->Pos.Z,
-	     simp->V[2]->Pos.X, simp->V[2]->Pos.Y, simp->V[2]->Pos.Z);
-    */
     prosca(n1, n2, &res);
     if(res < 0.0){
-      Msg(DEBUG, "Inverting orientation of elements in %s surface %d (res = %g)",
+      Msg(DEBUG, "Inverting triangles in %s surface %d (res = %g)",
 	  (s->Typ == MSH_SURF_PLAN) ? "Plane" : "NonPlane", s->Num, res);
-      Tree_Action(s->Simplexes, ActionInvertTriQua);
+      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);
+      Tree_Action(s->Quadrangles, ActionInvertQua);
     }
   }
 }
@@ -929,17 +940,11 @@ void Maillage_Surface(void *data, void *dum)
   if(MeshTransfiniteSurface(s) ||
      MeshEllipticSurface(s) ||
      MeshCylindricalSurface(s) ||
-     MeshParametricSurface(s)) {
-    Tree_Action(THEM->Points, PutVertex_OnSurf);
-    Tree_Action(s->Vertices, PutVertex_OnSurf);
-    Tree_Action(s->Vertices, Add_In_Mesh);
-  }
-  else if(Extrude_Mesh(s)) {
+     MeshParametricSurface(s) ||
+     Extrude_Mesh(s)) {
     Tree_Action(THEM->Points, PutVertex_OnSurf);
     Tree_Action(s->Vertices, PutVertex_OnSurf);
     Tree_Action(s->Vertices, Add_In_Mesh);
-    // FIXME: big hack
-    //ReOrientSurfaceMesh(s);
   }
   else{
     int TypSurface = s->Typ;

Mesh/3D_Extrude.cpp

-// $Id: 3D_Extrude.cpp,v 1.80 2004-06-08 02:10:32 geuzaine Exp $
+// $Id: 3D_Extrude.cpp,v 1.81 2004-06-23 03:57:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -1046,6 +1046,9 @@ 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,6 +476,9 @@ 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);
 void ActionLiss(void *data, void *dummy);

TODO

-$Id: TODO,v 1.49 2004-06-02 02:22:13 geuzaine Exp $
+$Id: TODO,v 1.50 2004-06-23 03:57:43 geuzaine Exp $
 
 add ternary operator and <,>,<=,>=,== tests in MathEval
 
@@ -34,6 +34,12 @@ interactive use?)
 
 ********************************************************************
 
+Find a better solution for the orientation of the surface mesh
+elements. The current "a posteriori re-orientation" solution is a
+real crime against good taste :-)
+
+********************************************************************
+
 Yves Krahenbuhl wrote:
 
 > Lors de la creation des elements du 2eme ordre, et selon la courbure
@@ -45,7 +51,7 @@ Yves Krahenbuhl wrote:
 
 ********************************************************************
 
-The "Symmetry" operation should be renamed "Reflection"
+The "Symmetry" operation should be renamed "Reflection" (?)
 
 ********************************************************************
 
@@ -55,7 +61,8 @@ Attractors in the 2D aniso algo are extremely buggy
 
 Memory leaks, memory leaks
 
-- start with mesh_domain() and the parser
+- start with mesh_domain() and the parser (update: parser should be
+  mostly OK now)
 
 - check all calls to Tree_Replace: we shouldn't use it with trees of
   pointers, since we loose the original pointer when we actually do a
@@ -90,6 +97,8 @@ etc.
 Include Tetgen (the same way we did it with Triangle--the license is
 the same).
 
+********************************************************************
+
 Include the 3D frontal algo from Netgen? (it's LGPL now...)
 
 ********************************************************************
@@ -162,68 +171,3 @@ Remarque finale: le calcul des valeurs et vecteurs propres est
 à partir des v.p.), et cela nous rappelle cette belle expression "No 
 free meal" !
 
-********************************************************************
-
-Stockage et acces aux normales en post-pro: generer les noeuds pour
-acces plus rapide?
-
-Idem pour le maillage.
-
-********************************************************************
-
-Orientation des surfaces
-
-mesh: il faut vraiment que je regarde la relation entre l'orientation
-des surfaces geometriques et les l'orientation des elements du
-maillage...
-
-Subject: Re: [Gmsh] surface normals: some point inwards, but most
-point in the negative direction
-        Date: Thu, 21 Nov 2002 05:59:56 -0500
-        From: "Matthias G. Imhof" <mgi@vt.edu>
- Organization: Virginia Tech: Seismic Reservoir Characterization
-Laboratory
-          To: Christophe Geuzaine <geuzaine@acm.caltech.edu>
-  References: 1 , 2
-
-Christophe,
-
-Thanks for a nice tool! I finally extended my boundary-element method
-to 3D and find your mesh generator perfect to make sure the boundary
-points are placed uniformly on a sphere instead of clustering up at
-the poles.
-
-My only gripe is with the visualization of surface normals in
-gmsh. While Physical Surface allows change of the normal directions in
-the mesh file, the normals are not changed in the visualization which
-makes identification of wrong polarity surfaces difficult.
-
-Matthias
-
-Christophe Geuzaine wrote:
-> Matthias G. Imhof wrote:
-> 
->> I defined a simple rectangular block and noticed that normals in the 
->> x-direction actually change signs depending on the surface to point 
->> into the block, while the other normals point in the negative y- or 
->> z-directions pointing into the block for some surfaces and pointing 
->> outwards for others.
->>
->> How can I change this model to have all normals point out of the block?
-> 
-> 
-> The only way is to define Physical Surfaces with the appropriate sign, 
-> e.g. "Physical Surface(1) = {23,1,20,3,2,-16};". This will reverse the 
-> orientation of the surface 16 in the output mesh file.
-> 
->>
->> A second question for the same problem. If I generate a quantity by 
->> duplicate, how can I lates access this quantity? I tried to assign it 
->> to, e.g., a plane surface but that did not work.
-> 
-> 
-> There is unfortunately no way to access that information at the moment 
-> without using the GUI...
-> 
-> Christophe
->