From 873ac3d272c18094949d1e07e4b1611c3331e256 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 25 May 2004 23:16:33 +0000
Subject: [PATCH] Added final bit for 2nd order elements: the middle face nodes
 for quadrangular faces (i.e., for quadrangles, hexahedra, prisms and
 pyramids).

---
 Common/Makefile        |  14 +-
 Fltk/Makefile          |  44 +--
 Geo/Makefile           |  30 +-
 Graphics/Makefile      |  73 ++---
 Graphics/Mesh.cpp      | 616 ++++++++++++++++++-----------------------
 Mesh/2D_Links.cpp      |   6 +-
 Mesh/2D_Mesh_Aniso.cpp |   6 +-
 Mesh/2D_Recombine.cpp  |  18 +-
 Mesh/3D_Coherence.cpp  |  20 +-
 Mesh/Create.cpp        |   6 +-
 Mesh/Edge.cpp          |  12 +-
 Mesh/Edge.h            |   5 +-
 Mesh/Face.cpp          | 287 +++++++++++++++++++
 Mesh/Face.h            |  65 +++++
 Mesh/Makefile          | 138 ++++-----
 Mesh/Mesh.h            |   9 +-
 Mesh/Print_Mesh.cpp    |  23 +-
 Mesh/Read_Mesh.cpp     |  26 +-
 Mesh/SecondOrder.cpp   | 174 ++++++++++--
 Mesh/Simplex.cpp       |  29 +-
 Mesh/Simplex.h         |   5 +-
 Mesh/SwapEdge.cpp      |   8 +-
 Parser/Gmsh.tab.cpp    | 301 ++++++++++----------
 Parser/Gmsh.y          |   3 +-
 Parser/Gmsh.yy.cpp     |   4 +-
 Parser/Makefile        |  12 +-
 Plugin/Makefile        |   4 +-
 doc/VERSIONS           |  27 +-
 doc/texinfo/gmsh.texi  | 138 ++++-----
 29 files changed, 1233 insertions(+), 870 deletions(-)
 create mode 100644 Mesh/Face.cpp
 create mode 100644 Mesh/Face.h

diff --git a/Common/Makefile b/Common/Makefile
index 9cc4e242b5..27d9a8b595 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.55 2004-04-24 03:51:59 geuzaine Exp $
+# $Id: Makefile,v 1.56 2004-05-25 23:16:25 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -62,7 +62,7 @@ depend:
 Context.o: Context.cpp Gmsh.h Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Numeric/Numeric.h ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Graphics/Draw.h ../Common/Views.h ../Common/ColorTable.h Context.h \
   Options.h DefaultOptions.h Trackball.h
@@ -74,15 +74,15 @@ Options.o: Options.cpp ../Plugin/PluginManager.h ../Plugin/Plugin.h \
   ../Common/ColorTable.h ../DataStr/List.h Gmsh.h ../DataStr/Malloc.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h GmshUI.h \
   ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Graphics/Draw.h Context.h \
-  ../Fltk/Solvers.h ../Fltk/GUI.h ../Fltk/Opengl_Window.h \
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Graphics/Draw.h \
+  Context.h ../Fltk/Solvers.h ../Fltk/GUI.h ../Fltk/Opengl_Window.h \
   ../Fltk/Colorbar_Window.h
 CommandLine.o: CommandLine.cpp Gmsh.h Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   GmshUI.h GmshVersion.h CommandLine.h ../Numeric/Numeric.h Context.h \
   Options.h ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   Views.h ColorTable.h ../Parser/OpenFile.h ../Parser/Parser.h
 Timer.o: Timer.cpp
@@ -92,7 +92,7 @@ ColorTable.o: ColorTable.cpp Gmsh.h Message.h ../DataStr/Malloc.h \
 Visibility.o: Visibility.cpp Gmsh.h Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Geo/Geo.h ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Parser/Parser.h Visibility.h
 Trackball.o: Trackball.cpp Trackball.h
diff --git a/Fltk/Makefile b/Fltk/Makefile
index 5b6164c438..cfd7bd82dd 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.52 2004-04-18 03:36:06 geuzaine Exp $
+# $Id: Makefile,v 1.53 2004-05-25 23:16:25 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -69,8 +69,8 @@ Main.o: Main.cpp ../Plugin/PluginManager.h ../Plugin/Plugin.h \
   ../DataStr/Malloc.h ../DataStr/Tree.h ../DataStr/avl.h \
   ../DataStr/Tools.h ../Common/GmshUI.h ../Common/GmshVersion.h \
   ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Graphics/Draw.h \
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Graphics/Draw.h \
   ../Common/Context.h ../Parser/Parser.h GUI.h Opengl_Window.h \
   Colorbar_Window.h ../Parser/OpenFile.h ../Common/CommandLine.h \
   ../Numeric/Numeric.h
@@ -85,17 +85,17 @@ GUI.o: GUI.cpp ../Plugin/PluginManager.h ../Plugin/Plugin.h \
   ../DataStr/Malloc.h ../DataStr/Tree.h ../DataStr/avl.h \
   ../DataStr/Tools.h ../Common/GmshUI.h ../Numeric/Numeric.h \
   ../Common/GmshVersion.h ../Common/Context.h ../Geo/Geo.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  ../Graphics/Draw.h GUI.h Opengl_Window.h Colorbar_Window.h Callbacks.h \
-  Bitmaps.h Win32Icon.h ../Parser/OpenFile.h ../Common/CommandLine.h \
-  Solvers.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h ../Graphics/Draw.h GUI.h Opengl_Window.h \
+  Colorbar_Window.h Callbacks.h Bitmaps.h Win32Icon.h \
+  ../Parser/OpenFile.h ../Common/CommandLine.h Solvers.h
 Callbacks.o: Callbacks.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Geo/ExtractContour.h \
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Geo/ExtractContour.h \
   ../Graphics/Draw.h ../Common/Views.h ../Common/ColorTable.h \
   ../Common/Timer.h ../Graphics/CreateFile.h ../Parser/OpenFile.h \
   ../Common/CommandLine.h ../Common/Context.h ../Common/Options.h GUI.h \
@@ -105,18 +105,19 @@ Opengl.o: Opengl.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h \
   ../Numeric/Numeric.h ../Common/Context.h ../Geo/Geo.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  ../Graphics/Draw.h ../Common/Views.h ../Common/ColorTable.h GUI.h \
-  Opengl_Window.h Colorbar_Window.h ../Graphics/gl2ps.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h ../Graphics/Draw.h ../Common/Views.h \
+  ../Common/ColorTable.h GUI.h Opengl_Window.h Colorbar_Window.h \
+  ../Graphics/gl2ps.h
 Opengl_Window.o: Opengl_Window.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h \
   ../Common/GmshUI.h ../Common/Context.h ../Geo/Geo.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  ../Graphics/Draw.h ../Common/Views.h ../Common/ColorTable.h GUI.h \
-  Opengl_Window.h Colorbar_Window.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h ../Graphics/Draw.h ../Common/Views.h \
+  ../Common/ColorTable.h GUI.h Opengl_Window.h Colorbar_Window.h
 Colorbar_Window.o: Colorbar_Window.cpp ../Common/Gmsh.h \
   ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
@@ -128,6 +129,7 @@ Solvers.o: Solvers.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../utils/solvers/GmshClient.h \
   GmshServer.h ../Parser/OpenFile.h Solvers.h ../Common/GmshUI.h GUI.h \
   Opengl_Window.h Colorbar_Window.h ../Common/ColorTable.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  ../Graphics/Draw.h ../Common/Views.h ../Common/Context.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h ../Graphics/Draw.h ../Common/Views.h \
+  ../Common/Context.h
diff --git a/Geo/Makefile b/Geo/Makefile
index 9604f4ab34..d8dd984abc 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.49 2004-05-25 04:10:03 geuzaine Exp $
+# $Id: Makefile,v 1.50 2004-05-25 23:16:26 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -61,7 +61,7 @@ depend:
 CAD.o: CAD.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Numeric/Numeric.h Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Mesh/Interpolation.h ../Mesh/Create.h CAD.h ../Common/Visibility.h \
   ../Common/Context.h
@@ -72,36 +72,38 @@ MinMax.o: MinMax.cpp ../Common/Gmsh.h ../Common/Message.h \
 ExtrudeParams.o: ExtrudeParams.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Geo.h CAD.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h
 Geo.o: Geo.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Numeric/Numeric.h Geo.h CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Parser/Parser.h ../Common/Context.h
 GeoUtils.o: GeoUtils.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Geo.h CAD.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h
 StepGeomDatabase.o: StepGeomDatabase.cpp ../Common/Gmsh.h \
   ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Numeric/Numeric.h Geo.h GeoUtils.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Mesh/Nurbs.h CAD.h StepGeomDatabase.h ../Mesh/Create.h \
   ../Common/Context.h
 ExtractContour.o: ExtractContour.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Geo.h GeoUtils.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  CAD.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h CAD.h
 Print_Geo.o: Print_Geo.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Geo.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  CAD.h ../Common/Context.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h CAD.h ../Common/Context.h
diff --git a/Graphics/Makefile b/Graphics/Makefile
index efc93e2aec..2dab2cbc3a 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.57 2004-05-17 17:42:24 geuzaine Exp $
+# $Id: Makefile,v 1.58 2004-05-25 23:16:26 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -72,43 +72,43 @@ depend:
 Draw.o: Draw.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Geo/Geo.h ../Geo/CAD.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  Draw.h ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
-  ../Geo/MinMax.h ../Numeric/Numeric.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h Draw.h ../Common/Views.h ../Common/ColorTable.h \
+  ../Common/Context.h ../Geo/MinMax.h ../Numeric/Numeric.h
 Mesh.o: Mesh.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Geo/Geo.h ../Geo/CAD.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  Draw.h ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
-  ../Geo/MinMax.h gl2ps.h ../Numeric/Numeric.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h Draw.h ../Common/Views.h ../Common/ColorTable.h \
+  ../Common/Context.h ../Geo/MinMax.h gl2ps.h ../Numeric/Numeric.h
 Geom.o: Geom.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Numeric/Numeric.h ../Geo/Geo.h ../Geo/CAD.h \
   ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h \
-  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
-  ../Mesh/Matrix.h ../Mesh/Utils.h Draw.h ../Common/Views.h \
-  ../Common/ColorTable.h ../Common/Context.h ../Mesh/Interpolation.h \
-  gl2ps.h
+  ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
+  ../Mesh/Metric.h ../Mesh/Matrix.h ../Mesh/Utils.h Draw.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
+  ../Mesh/Interpolation.h gl2ps.h
 Post.o: Post.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Numeric/Numeric.h ../Geo/Geo.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  Draw.h ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
-  gl2ps.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h Draw.h ../Common/Views.h ../Common/ColorTable.h \
+  ../Common/Context.h gl2ps.h
 PostElement.o: PostElement.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h ../Geo/Geo.h \
   ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h \
-  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
-  ../Mesh/Matrix.h Draw.h ../Common/Views.h ../Common/ColorTable.h Iso.h \
-  ../Common/Context.h ../Numeric/Numeric.h
+  ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
+  ../Mesh/Metric.h ../Mesh/Matrix.h Draw.h ../Common/Views.h \
+  ../Common/ColorTable.h Iso.h ../Common/Context.h ../Numeric/Numeric.h
 Iso.o: Iso.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   Draw.h ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
   ../Numeric/Numeric.h
@@ -116,9 +116,9 @@ Entity.o: Entity.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h \
   ../Numeric/Numeric.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Draw.h ../Common/Views.h \
-  ../Common/ColorTable.h ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Draw.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h
 ReadImg.o: ReadImg.cpp ReadImg.h ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h \
@@ -127,31 +127,32 @@ Scale.o: Scale.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h \
   ../Numeric/Numeric.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Draw.h ../Common/Views.h \
-  ../Common/ColorTable.h ../Common/Context.h gl2ps.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Draw.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h gl2ps.h
 Graph2D.o: Graph2D.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h \
   ../Common/Context.h ../Numeric/Numeric.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  Draw.h ../Common/Views.h ../Common/ColorTable.h gl2ps.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h Draw.h ../Common/Views.h ../Common/ColorTable.h \
+  gl2ps.h
 Axes.o: Axes.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../Common/GmshUI.h ../Numeric/Numeric.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   Draw.h ../Common/Views.h ../Common/ColorTable.h ../Common/Context.h \
   gl2ps.h
 CreateFile.o: CreateFile.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshUI.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
-  ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
-  ../Parser/OpenFile.h Draw.h ../Common/Views.h ../Common/ColorTable.h \
-  ../Common/Context.h ../Common/Options.h gl2ps.h gl2gif.h gl2jpeg.h \
-  gl2png.h gl2ppm.h gl2yuv.h
+  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h \
+  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
+  ../Mesh/Matrix.h ../Parser/OpenFile.h Draw.h ../Common/Views.h \
+  ../Common/ColorTable.h ../Common/Context.h ../Common/Options.h gl2ps.h \
+  gl2gif.h gl2jpeg.h gl2png.h gl2ppm.h gl2yuv.h
 gl2ps.o: gl2ps.cpp gl2ps.h
 gl2gif.o: gl2gif.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index e91d218bcf..898a113e01 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.83 2004-05-25 04:10:04 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.84 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -33,6 +33,11 @@
 extern Mesh *THEM;
 extern Context_T CTX;
 
+extern int edges_tetra[6][2];
+extern int edges_hexa[12][2];
+extern int edges_prism[9][2];
+extern int edges_pyramid[8][2];
+
 static DrawingColor theColor;
 static int thePhysical = 0;
 
@@ -407,24 +412,78 @@ void Draw_Mesh_Line(void *a, void *b)
   }
 }
 
-void glNormal3verts(Vertex *v0, Vertex *v1, Vertex *v2, double n[3])
+void _normal3points(double x0, double y0, double z0,
+		    double x1, double y1, double z1,
+		    double x2, double y2, double z2,
+		    double n[3])
 {
-  double x1x0, y1y0, z1z0, x2x0, y2y0, z2z0;
-  x1x0 = v1->Pos.X - v0->Pos.X;
-  y1y0 = v1->Pos.Y - v0->Pos.Y;
-  z1z0 = v1->Pos.Z - v0->Pos.Z;
-  x2x0 = v2->Pos.X - v0->Pos.X;
-  y2y0 = v2->Pos.Y - v0->Pos.Y;
-  z2z0 = v2->Pos.Z - v0->Pos.Z;
+  double x1x0 = x1 - x0;
+  double y1y0 = y1 - y0;
+  double z1z0 = z1 - z0;
+  double x2x0 = x2 - x0;
+  double y2y0 = y2 - y0;
+  double z2z0 = z2 - z0;
   n[0] = y1y0 * z2z0 - z1z0 * y2y0;
   n[1] = z1z0 * x2x0 - x1x0 * z2z0;
   n[2] = x1x0 * y2y0 - y1y0 * x2x0;
   glNormal3dv(n);
 }
 
+void _triFace(double x0, double y0, double z0,
+	      double x1, double y1, double z1,
+	      double x2, double y2, double z2)
+{
+  double n[3];
+  if(CTX.mesh.light) 
+    _normal3points(x0, y0, z0, x1, y1, z1, x2, y2, z2, n);
+  glVertex3d(x0, y0, z0);
+  glVertex3d(x1, y1, z1);
+  glVertex3d(x2, y2, z2);
+}
+
+void _triFace2(double *x, double *y, double *z,
+	       double *x2, double *y2, double *z2,
+	       int i0, int i1, int i2, 
+	       int j0, int j1, int j2)
+{
+  _triFace(x[i0], y[i0], z[i0], x2[j0], y2[j0], z2[j0], x2[j2], y2[j2], z2[j2]);
+  _triFace(x[i1], y[i1], z[i1], x2[j1], y2[j1], z2[j1], x2[j0], y2[j0], z2[j0]);
+  _triFace(x[i2], y[i2], z[i2], x2[j2], y2[j2], z2[j2], x2[j1], y2[j1], z2[j1]);
+  _triFace(x2[j0], y2[j0], z2[j0], x2[j1], y2[j1], z2[j1], x2[j2], y2[j2], z2[j2]);
+}
+
+void _quadFace(double *x, double *y, double *z,
+	       int i0, int i1, int i2, int i3)
+{
+  double n[3];
+  if(CTX.mesh.light) 
+    _normal3points(x[i0], y[i0], z[i0],
+		   x[i1], y[i1], z[i1],
+		   x[i2], y[i2], z[i2], n);
+  glVertex3d(x[i0], y[i0], z[i0]);
+  glVertex3d(x[i1], y[i1], z[i1]);
+  glVertex3d(x[i2], y[i2], z[i2]);
+  glVertex3d(x[i3], y[i3], z[i3]);
+}
+
+void _quadFace2(double *x, double *y, double *z,
+		double *x2, double *y2, double *z2,
+		int i0, int i1, int i2, int i3,
+		int j0, int j1, int j2, int j3, int j4)
+{
+  _triFace(x[i0], y[i0], z[i0], x2[j0], y2[j0], z2[j0], x2[j4], y2[j4], z2[j4]);
+  _triFace(x2[j0], y2[j0], z2[j0], x[i1], y[i1], z[i1], x2[j4], y2[j4], z2[j4]);
+  _triFace(x[i1], y[i1], z[i1], x2[j1], y2[j1], z2[j1], x2[j4], y2[j4], z2[j4]);
+  _triFace(x2[j1], y2[j1], z2[j1], x[i2], y[i2], z[i2], x2[j4], y2[j4], z2[j4]);
+  _triFace(x[i2], y[i2], z[i2], x2[j2], y2[j2], z2[j2], x2[j4], y2[j4], z2[j4]);
+  _triFace(x2[j2], y2[j2], z2[j2], x[i3], y[i3], z[i3], x2[j4], y2[j4], z2[j4]);
+  _triFace(x[i3], y[i3], z[i3], x2[j3], y2[j3], z2[j3], x2[j4], y2[j4], z2[j4]);
+  _triFace(x2[j3], y2[j3], z2[j3], x[i0], y[i0], z[i0], x2[j4], y2[j4], z2[j4]);
+}
+
 void Draw_Mesh_Triangle(void *a, void *b)
 {
-  double pX[6], pY[6], pZ[6];
+  double X[3], Y[3], Z[3], X2[3], Y2[3], Z2[3];
   double n[3];
   char Num[256];
 
@@ -447,19 +506,16 @@ void Draw_Mesh_Triangle(void *a, void *b)
   double Yc = (s->V[0]->Pos.Y + s->V[1]->Pos.Y + s->V[2]->Pos.Y) / 3.;
   double Zc = (s->V[0]->Pos.Z + s->V[1]->Pos.Z + s->V[2]->Pos.Z) / 3.;
 
-  int nn = s->VSUP ? 6 : 3;
-
-  int k = 0;
   for(int i = 0; i < 3; i++) {
-    pX[k] = Xc + CTX.mesh.explode * (s->V[i]->Pos.X - Xc);
-    pY[k] = Yc + CTX.mesh.explode * (s->V[i]->Pos.Y - Yc);
-    pZ[k] = Zc + CTX.mesh.explode * (s->V[i]->Pos.Z - Zc);
-    k++;
-    if(s->VSUP){
-      pX[k] = Xc + CTX.mesh.explode * (s->VSUP[i]->Pos.X - Xc);
-      pY[k] = Yc + CTX.mesh.explode * (s->VSUP[i]->Pos.Y - Yc);
-      pZ[k] = Zc + CTX.mesh.explode * (s->VSUP[i]->Pos.Z - Zc);
-      k++;
+    X[i] = Xc + CTX.mesh.explode * (s->V[i]->Pos.X - Xc);
+    Y[i] = Yc + CTX.mesh.explode * (s->V[i]->Pos.Y - Yc);
+    Z[i] = Zc + CTX.mesh.explode * (s->V[i]->Pos.Z - Zc);
+  }
+  if(s->VSUP){
+    for(int i = 0; i < 3; i++) {
+      X2[i] = Xc + CTX.mesh.explode * (s->VSUP[i]->Pos.X - Xc);
+      Y2[i] = Yc + CTX.mesh.explode * (s->VSUP[i]->Pos.Y - Yc);
+      Z2[i] = Zc + CTX.mesh.explode * (s->VSUP[i]->Pos.Z - Zc);
     }
   }
 
@@ -469,13 +525,10 @@ void Draw_Mesh_Triangle(void *a, void *b)
     glLineStipple(1, 0x0F0F);
     gl2psEnable(GL2PS_LINE_STIPPLE);
     glBegin(GL_LINES);
-    int incr = s->VSUP ? 2 : 1;
-    for(int i = 0; i < nn; i+=incr) {
+    for(int i = 0; i < 3; i++) {
       int j = i ? (i - 1) : 2;
       glVertex3d(Xc, Yc, Zc);
-      glVertex3d((pX[i] + pX[j]) / 2.,
-		 (pY[i] + pY[j]) / 2.,
-		 (pZ[i] + pZ[j]) / 2.);
+      glVertex3d((X[i] + X[j]) / 2., (Y[i] + Y[j]) / 2., (Z[i] + Z[j]) / 2.);
     }
     glEnd();
     glDisable(GL_LINE_STIPPLE);
@@ -483,7 +536,9 @@ void Draw_Mesh_Triangle(void *a, void *b)
   }
 
   if(CTX.mesh.normals || CTX.mesh.light)
-    glNormal3verts(s->V[0], s->V[1], s->V[2], n);
+    _normal3points(X[0], Y[0], Z[0], 
+		   X[1], Y[1], Z[1],
+		   X[2], Y[2], Z[2], n);
 
   if(CTX.mesh.surfaces_faces){
     glColor4ubv((GLubyte *) & CTX.color.mesh.line);
@@ -503,8 +558,10 @@ void Draw_Mesh_Triangle(void *a, void *b)
 
   if(CTX.mesh.surfaces_edges){
     glBegin(GL_LINE_LOOP);
-    for(int i = 0; i < nn; i++)
-      glVertex3d(pX[i], pY[i], pZ[i]);
+    for(int i = 0; i < 3; i++){
+      glVertex3d(X[i], Y[i], Z[i]);
+      if(s->VSUP) glVertex3d(X2[i], Y2[i], Z2[i]);
+    }
     glEnd();
   }
 
@@ -524,33 +581,18 @@ void Draw_Mesh_Triangle(void *a, void *b)
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     glEnable(GL_POLYGON_OFFSET_FILL);
 
-    if(!s->VSUP) { // first order elements
+    if(!s->VSUP) {
       glBegin(GL_TRIANGLES);
-      glVertex3d(pX[0], pY[0], pZ[0]);
-      glVertex3d(pX[1], pY[1], pZ[1]);
-      glVertex3d(pX[2], pY[2], pZ[2]);
+      glVertex3d(X[0], Y[0], Z[0]);
+      glVertex3d(X[1], Y[1], Z[1]);
+      glVertex3d(X[2], Y[2], Z[2]);
       glEnd();
     }
     else {
       glBegin(GL_TRIANGLES);
-      glNormal3verts(s->V[0], s->VSUP[0], s->VSUP[2], n);
-      glVertex3d(pX[0], pY[0], pZ[0]);
-      glVertex3d(pX[1], pY[1], pZ[1]);
-      glVertex3d(pX[5], pY[5], pZ[5]);
-      glNormal3verts(s->VSUP[0], s->V[1], s->VSUP[1], n);
-      glVertex3d(pX[1], pY[1], pZ[1]);
-      glVertex3d(pX[2], pY[2], pZ[2]);
-      glVertex3d(pX[3], pY[3], pZ[3]);
-      glNormal3verts(s->VSUP[1], s->V[2], s->VSUP[2], n);
-      glVertex3d(pX[3], pY[3], pZ[3]);
-      glVertex3d(pX[4], pY[4], pZ[4]);
-      glVertex3d(pX[5], pY[5], pZ[5]);
-      glNormal3verts(s->VSUP[0], s->VSUP[1], s->VSUP[2], n);
-      glVertex3d(pX[1], pY[1], pZ[1]);
-      glVertex3d(pX[3], pY[3], pZ[3]);
-      glVertex3d(pX[5], pY[5], pZ[5]);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 2, 0, 1, 2);
       glEnd();
-    }
+    }		    
     glDisable(GL_POLYGON_OFFSET_FILL);
     glDisable(GL_LIGHTING);
   }
@@ -577,7 +619,7 @@ void Draw_Mesh_Triangle(void *a, void *b)
 
 void Draw_Mesh_Quadrangle(void *a, void *b)
 {
-  double pX[8], pY[8], pZ[8];
+  double X[4], Y[4], Z[4], X2[5], Y2[5], Z2[5];
   double n[3];
   char Num[256];
 
@@ -600,19 +642,16 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
   double Zc = 0.25 * (q->V[0]->Pos.Z + q->V[1]->Pos.Z + 
 		      q->V[2]->Pos.Z + q->V[3]->Pos.Z);
 
-  int nn = q->VSUP ? 8 : 4;
-
-  int k = 0;
   for(int i = 0; i < 4; i++) {
-    pX[k] = Xc + CTX.mesh.explode * (q->V[i]->Pos.X - Xc);
-    pY[k] = Yc + CTX.mesh.explode * (q->V[i]->Pos.Y - Yc);
-    pZ[k] = Zc + CTX.mesh.explode * (q->V[i]->Pos.Z - Zc);
-    k++;
-    if(q->VSUP){
-      pX[k] = Xc + CTX.mesh.explode * (q->VSUP[i]->Pos.X - Xc);
-      pY[k] = Yc + CTX.mesh.explode * (q->VSUP[i]->Pos.Y - Yc);
-      pZ[k] = Zc + CTX.mesh.explode * (q->VSUP[i]->Pos.Z - Zc);
-      k++;
+    X[i] = Xc + CTX.mesh.explode * (q->V[i]->Pos.X - Xc);
+    Y[i] = Yc + CTX.mesh.explode * (q->V[i]->Pos.Y - Yc);
+    Z[i] = Zc + CTX.mesh.explode * (q->V[i]->Pos.Z - Zc);
+  }
+  if(q->VSUP){
+    for(int i = 0; i < 5; i++) {
+      X2[i] = Xc + CTX.mesh.explode * (q->VSUP[i]->Pos.X - Xc);
+      Y2[i] = Yc + CTX.mesh.explode * (q->VSUP[i]->Pos.Y - Yc);
+      Z2[i] = Zc + CTX.mesh.explode * (q->VSUP[i]->Pos.Z - Zc);
     }
   }
 
@@ -622,13 +661,10 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
     glLineStipple(1, 0x0F0F);
     gl2psEnable(GL2PS_LINE_STIPPLE);
     glBegin(GL_LINES);
-    int incr = q->VSUP ? 2 : 1;
-    for(int i = 0; i < nn; i+=incr) {
+    for(int i = 0; i < 4; i++) {
       int j = i ? (i - 1) : 3;
       glVertex3d(Xc, Yc, Zc);
-      glVertex3d((pX[i] + pX[j]) / 2.,
-		 (pY[i] + pY[j]) / 2.,
-		 (pZ[i] + pZ[j]) / 2.);
+      glVertex3d((X[i] + X[j]) / 2., (Y[i] + Y[j]) / 2., (Z[i] + Z[j]) / 2.);
     }
     glEnd();
     glDisable(GL_LINE_STIPPLE);
@@ -636,7 +672,9 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
   }
 
   if(CTX.mesh.normals || CTX.mesh.light)
-    glNormal3verts(q->V[0], q->V[1], q->V[2], n);
+    _normal3points(X[0], Y[0], Z[0], 
+		   X[1], Y[1], Z[1],
+		   X[2], Y[2], Z[2], n);
 
   if(CTX.mesh.surfaces_faces){
     glColor4ubv((GLubyte *) & CTX.color.mesh.line);
@@ -656,8 +694,11 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
 
   if(CTX.mesh.surfaces_edges){
     glBegin(GL_LINE_LOOP);
-    for(int i = 0; i < nn; i++)
-      glVertex3d(pX[i], pY[i], pZ[i]);
+    for(int i = 0; i < 4; i++){
+      glVertex3d(X[i], Y[i], Z[i]);
+      if(q->VSUP)
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+    }
     glEnd();
   }
 
@@ -676,19 +717,17 @@ void Draw_Mesh_Quadrangle(void *a, void *b)
 
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     glEnable(GL_POLYGON_OFFSET_FILL);
-    if(!q->VSUP) { // first order elements
+    if(!q->VSUP) {
       glBegin(GL_QUADS);
-      glVertex3d(pX[0], pY[0], pZ[0]);
-      glVertex3d(pX[1], pY[1], pZ[1]);
-      glVertex3d(pX[2], pY[2], pZ[2]);
-      glVertex3d(pX[3], pY[3], pZ[3]);
+      glVertex3d(X[0], Y[0], Z[0]);
+      glVertex3d(X[1], Y[1], Z[1]);
+      glVertex3d(X[2], Y[2], Z[2]);
+      glVertex3d(X[3], Y[3], Z[3]);
       glEnd();
     }
     else {
-      // FIXME: should subdivide...
-      glBegin(GL_POLYGON);
-      for(int i = 0; i < nn; i++)
-	glVertex3d(pX[i], pY[i], pZ[i]);
+      glBegin(GL_TRIANGLES);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 2, 3, 0, 1, 2, 3, 4);
       glEnd();
     }
     glDisable(GL_POLYGON_OFFSET_FILL);
@@ -795,37 +834,18 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
   }
 
   if(edges) {
-    if(!s->VSUP){
-      glBegin(GL_LINES);
-      glVertex3d(X[0], Y[0], Z[0]); glVertex3d(X[1], Y[1], Z[1]);
-      glVertex3d(X[1], Y[1], Z[1]); glVertex3d(X[2], Y[2], Z[2]);
-      glVertex3d(X[2], Y[2], Z[2]); glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X[2], Y[2], Z[2]); 
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X[1], Y[1], Z[1]);
-      glEnd();
-    }
-    else{
-      glBegin(GL_LINES);
-      glVertex3d(X[0], Y[0], Z[0]); glVertex3d(X2[0], Y2[0], Z2[0]); 
-      glVertex3d(X2[0], Y2[0], Z2[0]); glVertex3d(X[1], Y[1], Z[1]);
-
-      glVertex3d(X[1], Y[1], Z[1]); glVertex3d(X2[1], Y2[1], Z2[1]); 
-      glVertex3d(X2[1], Y2[1], Z2[1]); glVertex3d(X[2], Y[2], Z[2]);
-
-      glVertex3d(X[2], Y[2], Z[2]); glVertex3d(X2[2], Y2[2], Z2[2]); 
-      glVertex3d(X2[2], Y2[2], Z2[2]); glVertex3d(X[0], Y[0], Z[0]);
-
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X2[3], Y2[3], Z2[3]); 
-      glVertex3d(X2[3], Y2[3], Z2[3]); glVertex3d(X[0], Y[0], Z[0]);
-
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X2[4], Y2[4], Z2[4]); 
-      glVertex3d(X2[4], Y2[4], Z2[4]); glVertex3d(X[2], Y[2], Z[2]); 
-
-      glVertex3d(X[3], Y[3], Z[3]); glVertex3d(X2[5], Y2[5], Z2[5]); 
-      glVertex3d(X2[5], Y2[5], Z2[5]); glVertex3d(X[1], Y[1], Z[1]);
-      glEnd();
+    glBegin(GL_LINES);
+    for(int i = 0; i < 6; i++){
+      int j = edges_tetra[i][0];
+      int k = edges_tetra[i][1];
+      glVertex3d(X[j], Y[j], Z[j]);
+      if(s->VSUP){
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+      }
+      glVertex3d(X[k], Y[k], Z[k]);
     }
+    glEnd();
   }
 
   if(CTX.mesh.volumes_num) {
@@ -869,99 +889,22 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
     else
       glColor4ubv((GLubyte *) & CTX.color.mesh.tetrahedron);
 
-    double n[3];
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     if(CTX.mesh.surfaces_edges || edges) glEnable(GL_POLYGON_OFFSET_FILL);
     if(!s->VSUP){
       glBegin(GL_TRIANGLES);
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->V[2], s->V[1], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X[2], Y[2], Z[2]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->V[1], s->V[3], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      glVertex3d(X[3], Y[3], Z[3]);
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->V[3], s->V[2], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X[3], Y[3], Z[3]);
-      glVertex3d(X[2], Y[2], Z[2]);
-      if(CTX.mesh.light) glNormal3verts(s->V[3], s->V[1], s->V[2], n);
-      glVertex3d(X[3], Y[3], Z[3]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      glVertex3d(X[2], Y[2], Z[2]);
+      _triFace(X[0], Y[0], Z[0], X[2], Y[2], Z[2], X[1], Y[1], Z[1]);
+      _triFace(X[0], Y[0], Z[0], X[1], Y[1], Z[1], X[3], Y[3], Z[3]);
+      _triFace(X[0], Y[0], Z[0], X[3], Y[3], Z[3], X[2], Y[2], Z[2]);
+      _triFace(X[3], Y[3], Z[3], X[1], Y[1], Z[1], X[2], Y[2], Z[2]);
       glEnd();
     }
     else{
       glBegin(GL_TRIANGLES);
-      // face 1
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->VSUP[2], s->VSUP[0], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[0], s->VSUP[1], s->V[1], n);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[2], s->V[2], s->VSUP[1], n);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      glVertex3d(X[2], Y[2], Z[2]);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[0], s->VSUP[2], s->VSUP[1], n);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      // face 2
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->VSUP[0], s->VSUP[3], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[0], s->V[1], s->VSUP[5], n);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[3], s->VSUP[5], s->V[3], n);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
-      glVertex3d(X[3], Y[3], Z[3]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[0], s->VSUP[5], s->VSUP[3], n);
-      glVertex3d(X2[0], Y2[0], Z2[0]);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      // face 3
-      if(CTX.mesh.light) glNormal3verts(s->V[0], s->VSUP[3], s->VSUP[2], n);
-      glVertex3d(X[0], Y[0], Z[0]);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[3], s->V[3], s->VSUP[4], n);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      glVertex3d(X[3], Y[3], Z[3]);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[2], s->VSUP[4], s->V[2], n);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      glVertex3d(X[2], Y[2], Z[2]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[2], s->VSUP[3], s->VSUP[4], n);
-      glVertex3d(X2[2], Y2[2], Z2[2]);
-      glVertex3d(X2[3], Y2[3], Z2[3]);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      // face 4
-      if(CTX.mesh.light) glNormal3verts(s->V[3], s->VSUP[5], s->VSUP[4], n);
-      glVertex3d(X[3], Y[3], Z[3]);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[5], s->V[1], s->VSUP[1], n);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
-      glVertex3d(X[1], Y[1], Z[1]);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[4], s->VSUP[1], s->V[2], n);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      glVertex3d(X[2], Y[2], Z[2]);
-      if(CTX.mesh.light) glNormal3verts(s->VSUP[1], s->VSUP[4], s->VSUP[5], n);
-      glVertex3d(X2[1], Y2[1], Z2[1]);
-      glVertex3d(X2[4], Y2[4], Z2[4]);
-      glVertex3d(X2[5], Y2[5], Z2[5]);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 2, 1, 2, 1, 0);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 3, 0, 5, 3);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 3, 2, 3, 4, 2);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 3, 1, 2, 5, 1, 4);
       glEnd();
     }
     glDisable(GL_POLYGON_OFFSET_FILL);
@@ -971,12 +914,11 @@ void Draw_Mesh_Tetrahedron(void *a, void *b)
 
 void Draw_Mesh_Hexahedron(void *a, void *b)
 {
-  Hexahedron *h;
-  int i;
-  double Xc = 0.0, Yc = 0.0, Zc = 0.0, X[8], Y[8], Z[8];
+  double Xc = 0.0, Yc = 0.0, Zc = 0.0;
+  double X[8], Y[8], Z[8], X2[18], Y2[18], Z2[18];
   char Num[100];
 
-  h = *(Hexahedron **) a;
+  Hexahedron *h = *(Hexahedron **) a;
 
   if(!(h->Visible & VIS_MESH))
     return;
@@ -1002,7 +944,7 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
   if(intersectCutPlane(8, h->V, &edges, &faces) < 0)
     return;
 
-  for(i = 0; i < 8; i++) {
+  for(int i = 0; i < 8; i++) {
     Xc += h->V[i]->Pos.X;
     Yc += h->V[i]->Pos.Y;
     Zc += h->V[i]->Pos.Z;
@@ -1027,34 +969,31 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
       glColor4ubv((GLubyte *) & CTX.color.mesh.hexahedron);
   }
 
-  for(i = 0; i < 8; i++) {
+  for(int i = 0; i < 8; i++) {
     X[i] = Xc + CTX.mesh.explode * (h->V[i]->Pos.X - Xc);
     Y[i] = Yc + CTX.mesh.explode * (h->V[i]->Pos.Y - Yc);
     Z[i] = Zc + CTX.mesh.explode * (h->V[i]->Pos.Z - Zc);
   }
+  if(h->VSUP){
+    for(int i = 0; i < 18; i++) {
+      X2[i] = Xc + CTX.mesh.explode * (h->VSUP[i]->Pos.X - Xc);
+      Y2[i] = Yc + CTX.mesh.explode * (h->VSUP[i]->Pos.Y - Yc);
+      Z2[i] = Zc + CTX.mesh.explode * (h->VSUP[i]->Pos.Z - Zc);
+    }
+  }
 
   if(edges){
-    glBegin(GL_LINE_LOOP);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glEnd();
-    glBegin(GL_LINE_LOOP);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[6], Y[6], Z[6]);
-    glVertex3d(X[7], Y[7], Z[7]);
-    glEnd();
     glBegin(GL_LINES);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[6], Y[6], Z[6]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[7], Y[7], Z[7]);
+    for(int i = 0; i < 12; i++){
+      int j = edges_hexa[i][0];
+      int k = edges_hexa[i][1];
+      glVertex3d(X[j], Y[j], Z[j]);
+      if(h->VSUP){
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+      }
+      glVertex3d(X[k], Y[k], Z[k]);
+    }
     glEnd();
   }
 
@@ -1117,41 +1056,28 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
     else
       glColor4ubv((GLubyte *) & CTX.color.mesh.hexahedron);
 
-    double n[3];
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     if(CTX.mesh.surfaces_edges || edges) glEnable(GL_POLYGON_OFFSET_FILL);
-    glBegin(GL_QUADS);
-    if(CTX.mesh.light) glNormal3verts(h->V[0], h->V[2], h->V[1], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    if(CTX.mesh.light) glNormal3verts(h->V[4], h->V[5], h->V[6], n);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[6], Y[6], Z[6]);
-    glVertex3d(X[7], Y[7], Z[7]);
-    if(CTX.mesh.light) glNormal3verts(h->V[0], h->V[1], h->V[5], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    if(CTX.mesh.light) glNormal3verts(h->V[1], h->V[2], h->V[6], n);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[6], Y[6], Z[6]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    if(CTX.mesh.light) glNormal3verts(h->V[2], h->V[3], h->V[7], n);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[7], Y[7], Z[7]);
-    glVertex3d(X[6], Y[6], Z[6]);
-    if(CTX.mesh.light) glNormal3verts(h->V[0], h->V[4], h->V[7], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[7], Y[7], Z[7]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glEnd();
+    if(!h->VSUP){
+      glBegin(GL_QUADS);
+      _quadFace(X, Y, Z, 0, 3, 2, 1);
+      _quadFace(X, Y, Z, 4, 5, 6, 7);
+      _quadFace(X, Y, Z, 0, 1, 5, 4);
+      _quadFace(X, Y, Z, 1, 2, 6, 5);
+      _quadFace(X, Y, Z, 2, 3, 7, 6);
+      _quadFace(X, Y, Z, 0, 4, 7, 3);
+      glEnd();
+    }
+    else{
+      glBegin(GL_TRIANGLES);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 3, 2, 1, 1, 5, 3, 0, 12);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 4, 5, 6, 7, 8, 10, 11, 9, 17);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 5, 4, 0, 4, 8, 2, 13);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 1, 2, 6, 5, 3, 6, 10, 4, 15);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 2, 3, 7, 6, 5, 7, 11, 6, 16);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 4, 7, 3, 2, 9, 7, 1, 14);
+      glEnd();
+    }
     glDisable(GL_POLYGON_OFFSET_FILL);
     glDisable(GL_LIGHTING);
   }
@@ -1160,12 +1086,11 @@ void Draw_Mesh_Hexahedron(void *a, void *b)
 
 void Draw_Mesh_Prism(void *a, void *b)
 {
-  Prism *p;
-  int i;
-  double Xc = 0.0, Yc = 0.0, Zc = 0.0, X[6], Y[6], Z[6];
+  double Xc = 0.0, Yc = 0.0, Zc = 0.0;
+  double X[6], Y[6], Z[6], X2[12], Y2[12], Z2[12];
   char Num[100];
 
-  p = *(Prism **) a;
+  Prism *p = *(Prism **) a;
 
   if(!(p->Visible & VIS_MESH))
     return;
@@ -1191,7 +1116,7 @@ void Draw_Mesh_Prism(void *a, void *b)
   if(intersectCutPlane(6, p->V, &edges, &faces) < 0)
     return;
 
-  for(i = 0; i < 6; i++) {
+  for(int i = 0; i < 6; i++) {
     Xc += p->V[i]->Pos.X;
     Yc += p->V[i]->Pos.Y;
     Zc += p->V[i]->Pos.Z;
@@ -1216,30 +1141,31 @@ void Draw_Mesh_Prism(void *a, void *b)
       glColor4ubv((GLubyte *) & CTX.color.mesh.prism);
   }
 
-  for(i = 0; i < 6; i++) {
+  for(int i = 0; i < 6; i++) {
     X[i] = Xc + CTX.mesh.explode * (p->V[i]->Pos.X - Xc);
     Y[i] = Yc + CTX.mesh.explode * (p->V[i]->Pos.Y - Yc);
     Z[i] = Zc + CTX.mesh.explode * (p->V[i]->Pos.Z - Zc);
   }
+  if(p->VSUP){
+    for(int i = 0; i < 12; i++) {
+      X2[i] = Xc + CTX.mesh.explode * (p->VSUP[i]->Pos.X - Xc);
+      Y2[i] = Yc + CTX.mesh.explode * (p->VSUP[i]->Pos.Y - Yc);
+      Z2[i] = Zc + CTX.mesh.explode * (p->VSUP[i]->Pos.Z - Zc);
+    }
+  }
 
   if(edges){
-    glBegin(GL_LINE_LOOP);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glEnd();
-    glBegin(GL_LINE_LOOP);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glEnd();
     glBegin(GL_LINES);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[5], Y[5], Z[5]);
+    for(int i = 0; i < 9; i++){
+      int j = edges_prism[i][0];
+      int k = edges_prism[i][1];
+      glVertex3d(X[j], Y[j], Z[j]);
+      if(p->VSUP){
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+      }
+      glVertex3d(X[k], Y[k], Z[k]);
+    }
     glEnd();
   }
 
@@ -1297,36 +1223,28 @@ void Draw_Mesh_Prism(void *a, void *b)
     else
       glColor4ubv((GLubyte *) & CTX.color.mesh.prism);
 
-    double n[3];
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     if(CTX.mesh.surfaces_edges || edges) glEnable(GL_POLYGON_OFFSET_FILL);
-    glBegin(GL_TRIANGLES);
-    if(CTX.mesh.light) glNormal3verts(p->V[0], p->V[2], p->V[1], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    if(CTX.mesh.light) glNormal3verts(p->V[3], p->V[4], p->V[5], n);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glEnd();
-    glBegin(GL_QUADS);
-    if(CTX.mesh.light) glNormal3verts(p->V[0], p->V[1], p->V[4], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    if(CTX.mesh.light) glNormal3verts(p->V[1], p->V[2], p->V[5], n);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    if(CTX.mesh.light) glNormal3verts(p->V[0], p->V[3], p->V[5], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[5], Y[5], Z[5]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glEnd();
+    if(!p->VSUP){
+      glBegin(GL_TRIANGLES);
+      _triFace(X[0], Y[0], Z[0], X[2], Y[2], Z[2], X[1], Y[1], Z[1]);
+      _triFace(X[3], Y[3], Z[3], X[4], Y[4], Z[4], X[5], Y[5], Z[5]);
+      glEnd();
+      glBegin(GL_QUADS);
+      _quadFace(X, Y, Z, 0, 1, 4, 3);
+      _quadFace(X, Y, Z, 1, 2, 5, 4);
+      _quadFace(X, Y, Z, 0, 3, 5, 2);
+      glEnd();
+    }
+    else{
+      glBegin(GL_TRIANGLES);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 2, 1, 1, 3, 0);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 3, 4, 5, 6, 8, 7);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 4, 3, 0, 4, 6, 2, 9);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 1, 2, 5, 4, 3, 5, 8, 4, 11);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 3, 5, 2, 2, 7, 5, 1, 10);
+      glEnd();
+    }
     glDisable(GL_POLYGON_OFFSET_FILL);
     glDisable(GL_LIGHTING);
   }
@@ -1334,12 +1252,11 @@ void Draw_Mesh_Prism(void *a, void *b)
 
 void Draw_Mesh_Pyramid(void *a, void *b)
 {
-  Pyramid *p;
-  int i;
-  double Xc = 0.0, Yc = 0.0, Zc = 0.0, X[5], Y[5], Z[5];
+  double Xc = 0.0, Yc = 0.0, Zc = 0.0;
+  double X[5], Y[5], Z[5], X2[9], Y2[9], Z2[9];
   char Num[100];
 
-  p = *(Pyramid **) a;
+  Pyramid *p = *(Pyramid **) a;
 
   if(!(p->Visible & VIS_MESH))
     return;
@@ -1365,7 +1282,7 @@ void Draw_Mesh_Pyramid(void *a, void *b)
   if(intersectCutPlane(5, p->V, &edges, &faces) < 0)
     return;
 
-  for(i = 0; i < 5; i++) {
+  for(int i = 0; i < 5; i++) {
     Xc += p->V[i]->Pos.X;
     Yc += p->V[i]->Pos.Y;
     Zc += p->V[i]->Pos.Z;
@@ -1390,28 +1307,31 @@ void Draw_Mesh_Pyramid(void *a, void *b)
       glColor4ubv((GLubyte *) & CTX.color.mesh.pyramid);
   }
 
-  for(i = 0; i < 5; i++) {
+  for(int i = 0; i < 5; i++) {
     X[i] = Xc + CTX.mesh.explode * (p->V[i]->Pos.X - Xc);
     Y[i] = Yc + CTX.mesh.explode * (p->V[i]->Pos.Y - Yc);
     Z[i] = Zc + CTX.mesh.explode * (p->V[i]->Pos.Z - Zc);
   }
+  if(p->VSUP){
+    for(int i = 0; i < 9; i++) {
+      X2[i] = Xc + CTX.mesh.explode * (p->VSUP[i]->Pos.X - Xc);
+      Y2[i] = Yc + CTX.mesh.explode * (p->VSUP[i]->Pos.Y - Yc);
+      Z2[i] = Zc + CTX.mesh.explode * (p->VSUP[i]->Pos.Z - Zc);
+    }
+  }
 
   if(edges){
-    glBegin(GL_LINE_LOOP);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glEnd();
     glBegin(GL_LINES);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[4], Y[4], Z[4]);
+    for(int i = 0; i < 8; i++){
+      int j = edges_pyramid[i][0];
+      int k = edges_pyramid[i][1];
+      glVertex3d(X[j], Y[j], Z[j]);
+      if(p->VSUP){
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+	glVertex3d(X2[i], Y2[i], Z2[i]);
+      }
+      glVertex3d(X[k], Y[k], Z[k]);
+    }
     glEnd();
   }
 
@@ -1433,34 +1353,28 @@ void Draw_Mesh_Pyramid(void *a, void *b)
     else
       glColor4ubv((GLubyte *) & CTX.color.mesh.pyramid);
 
-    double n[3];
     if(CTX.mesh.light) glEnable(GL_LIGHTING);
     if(CTX.mesh.surfaces_edges || edges) glEnable(GL_POLYGON_OFFSET_FILL);
-    glBegin(GL_QUADS);
-    if(CTX.mesh.light) glNormal3verts(p->V[0], p->V[3], p->V[2], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glEnd();
-    glBegin(GL_TRIANGLES);
-    if(CTX.mesh.light) glNormal3verts(p->V[1], p->V[2], p->V[4], n);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    if(CTX.mesh.light) glNormal3verts(p->V[2], p->V[3], p->V[4], n);
-    glVertex3d(X[2], Y[2], Z[2]);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    if(CTX.mesh.light) glNormal3verts(p->V[3], p->V[0], p->V[4], n);
-    glVertex3d(X[3], Y[3], Z[3]);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    if(CTX.mesh.light) glNormal3verts(p->V[0], p->V[1], p->V[4], n);
-    glVertex3d(X[0], Y[0], Z[0]);
-    glVertex3d(X[1], Y[1], Z[1]);
-    glVertex3d(X[4], Y[4], Z[4]);
-    glEnd();
+    if(!p->VSUP){
+      glBegin(GL_QUADS);
+      _quadFace(X, Y, Z, 0, 3, 2, 1);
+      glEnd();
+      glBegin(GL_TRIANGLES);
+      _triFace(X[1], Y[1], Z[1], X[2], Y[2], Z[2], X[4], Y[4], Z[4]);
+      _triFace(X[2], Y[2], Z[2], X[3], Y[3], Z[3], X[4], Y[4], Z[4]);
+      _triFace(X[3], Y[3], Z[3], X[0], Y[0], Z[0], X[4], Y[4], Z[4]);
+      _triFace(X[0], Y[0], Z[0], X[1], Y[1], Z[1], X[4], Y[4], Z[4]);
+      glEnd();
+    }
+    else{
+      glBegin(GL_TRIANGLES);
+      _quadFace2(X, Y, Z, X2, Y2, Z2, 0, 3, 2, 1, 1, 5, 3, 0, 8);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 1, 2, 4, 3, 6, 4);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 2, 3, 4, 5, 7, 6);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 3, 0, 4, 1, 2, 7);
+      _triFace2(X, Y, Z, X2, Y2, Z2, 0, 1, 4, 0, 4, 2);
+      glEnd();
+    }
     glDisable(GL_POLYGON_OFFSET_FILL);
     glDisable(GL_LIGHTING);
   }
diff --git a/Mesh/2D_Links.cpp b/Mesh/2D_Links.cpp
index 703a91a162..1214da9a3f 100644
--- a/Mesh/2D_Links.cpp
+++ b/Mesh/2D_Links.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Links.cpp,v 1.18 2004-03-03 22:26:33 geuzaine Exp $
+// $Id: 2D_Links.cpp,v 1.19 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -139,7 +139,7 @@ int Conversion(DocPeek doc)
 }
 
 
-int compareEdge(const void *e1, const void *e2)
+int compareedge(const void *e1, const void *e2)
 {
   PointNumero tmp;
 
@@ -207,7 +207,7 @@ int CreateLinks(List_T * ListDelaunay, int NumDelaunay,
 
   }
 
-  qsort(ListEdges, 3 * NumDelaunay, sizeof(edge), compareEdge);
+  qsort(ListEdges, 3 * NumDelaunay, sizeof(edge), compareedge);
 
   i = 0;
 
diff --git a/Mesh/2D_Mesh_Aniso.cpp b/Mesh/2D_Mesh_Aniso.cpp
index 76ce117b87..50ce7e9548 100644
--- a/Mesh/2D_Mesh_Aniso.cpp
+++ b/Mesh/2D_Mesh_Aniso.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh_Aniso.cpp,v 1.39 2004-05-25 04:10:04 geuzaine Exp $
+// $Id: 2D_Mesh_Aniso.cpp,v 1.40 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -291,7 +291,7 @@ int Intersect_Edges_2d(Edge * a, Edge * b)
   return 0;
 }
 
-int compareedgeptr(const void *a, const void *b)
+int compareEdgePtr(const void *a, const void *b)
 {
   int i1, i2, j1, j2;
   Edge *q, *w;
@@ -354,7 +354,7 @@ void Recover_Edge(Surface * s, Edge * e, EdgesContainer & Edges)
     List_Read(coquille, i, &me);
     v[0] = me->V[0];
     v[1] = me->V[1];
-    List_Suppress(coquille, &me, compareedgeptr);
+    List_Suppress(coquille, &me, compareEdgePtr);
     Edges.SwapEdge(v);
     if(Edges.Search(e->V[0], e->V[1]))
       break;
diff --git a/Mesh/2D_Recombine.cpp b/Mesh/2D_Recombine.cpp
index 46587078c5..28557a7406 100644
--- a/Mesh/2D_Recombine.cpp
+++ b/Mesh/2D_Recombine.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Recombine.cpp,v 1.19 2004-05-25 04:10:04 geuzaine Exp $
+// $Id: 2D_Recombine.cpp,v 1.20 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -34,8 +34,8 @@ static List_T *SimplexesToRemove;
 static double ALPHA;
 static int RECNUM;
 
-// Warning: these routines temporarily leave quads in the simplex
-// tree!
+// Note: these routines temporarily leave quads in the simplex tree,
+// and only remove them at the end...
 
 void addTriangles(void *a, void *b)
 {
@@ -88,7 +88,6 @@ void recombineFace(void *a, void *b)
       return;
     Tree_Add(RecSimplex, &s1);
     Tree_Suppress(TREEELM, &s1);
-    List_Add(SimplexesToRemove, &s1);
     s2->V[3] = ed->V[0];
     s2->V[2] = ed->O[0];
     s2->V[1] = ed->V[1];
@@ -118,10 +117,9 @@ int Recombine(Tree_T * Vertices, Tree_T * Simplexes, Tree_T * Quadrangles,
 
     // Initialization
     RECNUM = 0;
-    RecEdges = Tree_Create(sizeof(Edge), compareedge_angle);
+    RecEdges = Tree_Create(sizeof(Edge), compareEdgeAngle);
     RecSimplex = Tree_Create(sizeof(Simplex *), compareSimplex);
     Triangles = Tree_Create(sizeof(Simplex *), compareSimplex);
-    SimplexesToRemove = List_Create(100, 100, sizeof(Simplex*));
 
     // Recombination
     Tree_Action(Simplexes, addTriangles);
@@ -142,16 +140,10 @@ int Recombine(Tree_T * Vertices, Tree_T * Simplexes, Tree_T * Quadrangles,
 
     // Destruction
     Tree_Delete(RecEdges);
+    Tree_Action(RecSimplex, Free_Simplex);
     Tree_Delete(RecSimplex);
     Tree_Delete(Triangles);
 
-    for(int i = 0; i < List_Nbr(SimplexesToRemove); i++){
-      Simplex *s;
-      List_Read(SimplexesToRemove, i, &s); 
-      Free_Simplex(&s, NULL);
-    }
-    List_Delete(SimplexesToRemove);
-
     ntot += RECNUM;
     if(!RECNUM)
       break;
diff --git a/Mesh/3D_Coherence.cpp b/Mesh/3D_Coherence.cpp
index 3a1585402c..70201a7cb0 100644
--- a/Mesh/3D_Coherence.cpp
+++ b/Mesh/3D_Coherence.cpp
@@ -1,4 +1,4 @@
-// $Id: 3D_Coherence.cpp,v 1.33 2004-04-19 00:18:07 geuzaine Exp $
+// $Id: 3D_Coherence.cpp,v 1.34 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -97,7 +97,7 @@ void find_quads(void *a, void *b)
   Simplex *s1, *s2;
   q = (Edge *) a;
 
-  if(!List_Search(Missing, q, compareedge))
+  if(!List_Search(Missing, q, compareEdge))
     return;
 
   if(List_Nbr(q->Simplexes) != 2)
@@ -174,7 +174,7 @@ void swap_quads(void *a, void *b)
   qsort(s1->F[0].V, 3, sizeof(Vertex *), compareVertex);
   qsort(s2->F[0].V, 3, sizeof(Vertex *), compareVertex);
 
-  List_Suppress(Missing, q, compareedge);
+  List_Suppress(Missing, q, compareEdge);
   q->V[0] = q->O[0];
   q->V[1] = q->O[1];
 
@@ -184,7 +184,7 @@ int create_Quads(Volume * V)
 {
   int i, n;
   Surface *S;
-  swaps = Tree_Create(sizeof(Edge), compareedge);
+  swaps = Tree_Create(sizeof(Edge), compareEdge);
   touchedvertex = Tree_Create(sizeof(Vertex *), compareVertex);
   for(i = 0; i < List_Nbr(V->Surfaces); i++) {
     List_Read(V->Surfaces, i, &S);
@@ -279,7 +279,7 @@ void create_Edges(Volume * V)
     Tree_Delete(V->Edges);
   }
 
-  V->Edges = Tree_Create(sizeof(Edge), compareedge);
+  V->Edges = Tree_Create(sizeof(Edge), compareEdge);
   EdgesTree = V->Edges;
 
   Tree_Action(V->Simplexes, create_Edge);
@@ -291,7 +291,7 @@ void create_Edges(Volume * V)
       //Tree_Action(S->Edges,Free_Edge);
       Tree_Delete(S->Edges);
     }
-    S->Edges = Tree_Create(sizeof(Edge), compareedge);
+    S->Edges = Tree_Create(sizeof(Edge), compareEdge);
     EdgesTree = S->Edges;
     Tree_Action(S->Simplexes, create_Edge);
   }
@@ -418,7 +418,7 @@ int compxNewv(const void *a, const void *b)
   if(q->ef != w->ef)
     return (q->ef - w->ef);
   if(q->ef == 1)
-    return compareedge(&q->e, &w->e);
+    return compareEdge(&q->e, &w->e);
   if(q->ef == 2)
     return compareFace(q->f, w->f);
   return 1;
@@ -1246,13 +1246,13 @@ int Coherence(Volume * v, Mesh * m)
         simp->F[0].V[1]->Num, simp->F[0].V[2]->Num);
     E.V[0] = simp->F[0].V[0];
     E.V[1] = simp->F[0].V[1];
-    pE1 = (Edge *) List_PQuery(Missing, &E, compareedge);
+    pE1 = (Edge *) List_PQuery(Missing, &E, compareEdge);
     E.V[0] = simp->F[0].V[1];
     E.V[1] = simp->F[0].V[2];
-    pE2 = (Edge *) List_PQuery(Missing, &E, compareedge);
+    pE2 = (Edge *) List_PQuery(Missing, &E, compareEdge);
     E.V[0] = simp->F[0].V[2];
     E.V[1] = simp->F[0].V[0];
-    pE3 = (Edge *) List_PQuery(Missing, &E, compareedge);
+    pE3 = (Edge *) List_PQuery(Missing, &E, compareEdge);
 
     /* On verifie si c'est simple c a d si les tetraedres
        couvrent entierement la face */
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index 8f0390b745..9f02e1281a 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.55 2004-05-25 04:10:04 geuzaine Exp $
+// $Id: Create.cpp,v 1.56 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -41,7 +41,7 @@ int compareNXE(const void *a, const void *b)
 
   q = (NXE *) a;
   w = (NXE *) b;
-  return (compareVertex(&q->v, &w->v));
+  return compareVertex(&q->v, &w->v);
 }
 
 int compareFxE(const void *a, const void *b)
@@ -50,7 +50,7 @@ int compareFxE(const void *a, const void *b)
 
   q = (FxE *) a;
   w = (FxE *) b;
-  return (compareFace(&q->Sorted, &w->Sorted));
+  return compareFace(&q->Sorted, &w->Sorted);
 }
 
 int compareSurfaceLoop(const void *a, const void *b)
diff --git a/Mesh/Edge.cpp b/Mesh/Edge.cpp
index 2c8d22d270..d09740dec8 100644
--- a/Mesh/Edge.cpp
+++ b/Mesh/Edge.cpp
@@ -1,4 +1,4 @@
-// $Id: Edge.cpp,v 1.16 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Edge.cpp,v 1.17 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -81,7 +81,7 @@ int edges_pyramid[8][2] = {
 
 int edges_non[3] = { 2, 0, 1 };
 
-int compareedge(const void *a, const void *b)
+int compareEdge(const void *a, const void *b)
 {
   int i1, i2, j1, j2;
   Edge *q, *w;
@@ -104,7 +104,7 @@ int compareedge(const void *a, const void *b)
   return 0;
 }
 
-int compareedge_angle(const void *a, const void *b)
+int compareEdgeAngle(const void *a, const void *b)
 {
   Edge *q, *w;
 
@@ -236,13 +236,13 @@ void EdgesContainer::AddEdges(Pyramid * p)
 
 EdgesContainer::EdgesContainer(Tree_T * Simplexes)
 {
-  AllEdges = Tree_Create(sizeof(Edge), compareedge);
+  AllEdges = Tree_Create(sizeof(Edge), compareEdge);
   AddSimplexTree(Simplexes);
 }
 
 EdgesContainer::EdgesContainer(List_T * Surfaces)
 {
-  AllEdges = Tree_Create(sizeof(Edge), compareedge);
+  AllEdges = Tree_Create(sizeof(Edge), compareEdge);
   Surface *s;
   for(int i = 0; i < List_Nbr(Surfaces); i++) {
     List_Read(Surfaces, i, &s);
@@ -252,7 +252,7 @@ EdgesContainer::EdgesContainer(List_T * Surfaces)
 
 EdgesContainer::EdgesContainer()
 {
-  AllEdges = Tree_Create(sizeof(Edge), compareedge);
+  AllEdges = Tree_Create(sizeof(Edge), compareEdge);
 }
 
 void EdgesContainer::AddSimplexTree(Tree_T * Simplexes)
diff --git a/Mesh/Edge.h b/Mesh/Edge.h
index a71167ce8d..75241708cc 100644
--- a/Mesh/Edge.h
+++ b/Mesh/Edge.h
@@ -76,8 +76,9 @@ class EdgesContainer
     void Print();
 };
 
-int compareedge(const void *a, const void *b);
-int compareedge_angle(const void *a, const void *b);
 void Free_Edge(void *a, void *b);
 
+int compareEdge(const void *a, const void *b);
+int compareEdgeAngle(const void *a, const void *b);
+
 #endif
diff --git a/Mesh/Face.cpp b/Mesh/Face.cpp
new file mode 100644
index 0000000000..0fbc29177e
--- /dev/null
+++ b/Mesh/Face.cpp
@@ -0,0 +1,287 @@
+// $Id: Face.cpp,v 1.1 2004-05-25 23:16:26 geuzaine Exp $
+//
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Gmsh.h"
+#include "Numeric.h"
+#include "Mesh.h"
+
+// Simplex faces
+
+extern int FACE_DIMENSION;
+
+int compareFace(const void *a, const void *b)
+{
+  Face *q = (Face *) a;
+  Face *w = (Face *) b;
+
+  if(!q->V[0] || !w->V[0])
+    Msg(FATAL, "Bad face (are you trying to generate hexahedra with a Delaunay?!)");
+
+  if(q->V[0]->Num > w->V[0]->Num)
+    return (1);
+  if(q->V[0]->Num < w->V[0]->Num)
+    return (-1);
+
+  if(q->V[1]->Num > w->V[1]->Num)
+    return (1);
+  if(q->V[1]->Num < w->V[1]->Num)
+    return (-1);
+
+  if(FACE_DIMENSION == 1 || !q->V[2] || !w->V[2])
+    return 0;
+
+  if(q->V[2]->Num > w->V[2]->Num)
+    return (1);
+  if(q->V[2]->Num < w->V[2]->Num)
+    return (-1);
+  return (0);
+}
+
+// Quad faces
+
+int quadfaces_hexa[6][4] = {
+  {0, 1, 2, 3},
+  {0, 1, 4, 5},
+  {0, 3, 4, 7},
+  {1, 2, 5, 6},
+  {2, 3, 6, 7},
+  {4, 5, 6, 7}
+};
+
+int quadfaces_prism[3][4] = {
+  {0, 1, 3, 4},
+  {0, 2, 3, 5},
+  {1, 2, 4, 5}
+};
+
+QuadFace::QuadFace()
+{
+  V[0] = V[1] = V[2] = V[3] = NULL;
+  newv = NULL;
+  quad = NULL;
+  hexa[0] = hexa[1] = NULL;
+  prism[0] = prism[1] = NULL;
+  pyramid[0] = pyramid[1] = NULL;
+}
+
+QuadFacesContainer::QuadFacesContainer()
+{
+  AllFaces = Tree_Create(sizeof(QuadFace), compareQuadFace);
+}
+
+QuadFacesContainer::~QuadFacesContainer()
+{
+  Tree_Delete(AllFaces);
+}
+
+void QuadFacesContainer::AddQuadrangleTree(Tree_T *Quadrangles)
+{
+  Quadrangle *q;
+  List_T *temp = Tree2List(Quadrangles);
+  for(int i = 0; i < List_Nbr(temp); i++) {
+    List_Read(temp, i, &q);
+    AddQuadFaces(q);
+  }
+  List_Delete(temp);
+}
+
+void QuadFacesContainer::AddHexahedronTree(Tree_T *Hexahedra)
+{
+  Hexahedron *h;
+  List_T *temp = Tree2List(Hexahedra);
+  for(int i = 0; i < List_Nbr(temp); i++) {
+    List_Read(temp, i, &h);
+    AddQuadFaces(h);
+  }
+  List_Delete(temp);
+}
+
+void QuadFacesContainer::AddPrismTree(Tree_T *Prisms)
+{
+  Prism *p;
+  List_T *temp = Tree2List(Prisms);
+  for(int i = 0; i < List_Nbr(temp); i++) {
+    List_Read(temp, i, &p);
+    AddQuadFaces(p);
+  }
+  List_Delete(temp);
+}
+
+void QuadFacesContainer::AddPyramidTree(Tree_T *Pyramids)
+{
+  Pyramid *p;
+  List_T *temp = Tree2List(Pyramids);
+  for(int i = 0; i < List_Nbr(temp); i++) {
+    List_Read(temp, i, &p);
+    AddQuadFaces(p);
+  }
+  List_Delete(temp);
+}
+
+void QuadFacesContainer::AddQuadFaces(Quadrangle *q)
+{
+  QuadFace F, *pF;
+
+  F.V[0] = q->V[0];
+  F.V[1] = q->V[1];
+  F.V[2] = q->V[2];
+  F.V[3] = q->V[3];
+  qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+
+  if((pF = (QuadFace *) Tree_PQuery(AllFaces, &F))) {
+    if(!pF->quad)
+      pF->quad = q;
+  }
+  else {
+    F.quad = q;
+    F.newv = NULL;
+    Tree_Add(AllFaces, &F);
+  }
+}
+
+void QuadFacesContainer::AddQuadFaces(Hexahedron *h)
+{
+  QuadFace F, *pF;
+
+  for(int i = 0; i < 6; i++) {
+    F.V[0] = h->V[quadfaces_hexa[i][0]];
+    F.V[1] = h->V[quadfaces_hexa[i][1]];
+    F.V[2] = h->V[quadfaces_hexa[i][2]];
+    F.V[3] = h->V[quadfaces_hexa[i][3]];
+    qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+
+    if((pF = (QuadFace *) Tree_PQuery(AllFaces, &F))) {
+      if(!pF->hexa[0])
+	pF->hexa[0] = h;
+      else if(!pF->hexa[1])
+	pF->hexa[1] = h;
+    }
+    else {
+      F.hexa[0] = h;
+      F.newv = NULL;
+      Tree_Add(AllFaces, &F);
+    }
+  }
+}
+
+void QuadFacesContainer::AddQuadFaces(Prism *p)
+{
+  QuadFace F, *pF;
+
+  for(int i = 0; i < 3; i++) {
+    F.V[0] = p->V[quadfaces_prism[i][0]];
+    F.V[1] = p->V[quadfaces_prism[i][1]];
+    F.V[2] = p->V[quadfaces_prism[i][2]];
+    F.V[3] = p->V[quadfaces_prism[i][3]];
+    qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+
+    if((pF = (QuadFace *) Tree_PQuery(AllFaces, &F))) {
+      if(!pF->prism[0])
+	pF->prism[0] = p;
+      else if(!pF->prism[1])
+	pF->prism[1] = p;
+    }
+    else {
+      F.prism[0] = p;
+      F.newv = NULL;
+      Tree_Add(AllFaces, &F);
+    }
+  }
+}
+
+void QuadFacesContainer::AddQuadFaces(Pyramid *p)
+{
+  QuadFace F, *pF;
+
+  F.V[0] = p->V[0];
+  F.V[1] = p->V[1];
+  F.V[2] = p->V[2];
+  F.V[3] = p->V[3];
+  qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+
+  if((pF = (QuadFace *) Tree_PQuery(AllFaces, &F))) {
+    if(!pF->pyramid[0])
+      pF->pyramid[0] = p;
+    else if(!pF->pyramid[1])
+      pF->pyramid[1] = p;
+  }
+  else {
+    F.pyramid[0] = p;
+    F.newv = NULL;
+    Tree_Add(AllFaces, &F);
+  }
+}
+
+void QuadFacesContainer::Print()
+{
+  List_T *temp = Tree2List(AllFaces);
+  printf("Print QuadFaces START\n");
+  for(int i = 0; i < List_Nbr(temp); i++) {
+    QuadFace *f = (QuadFace*)List_Pointer(temp, i);
+    printf("quadface %d %d %d %d", 
+	   f->V[0]->Num, f->V[1]->Num, f->V[2]->Num, f->V[3]->Num);
+    if(f->quad)
+      printf(", in quad %d", f->quad->Num);
+    if(f->hexa[0])
+      printf(", in hexa %d", f->hexa[0]->Num);
+    if(f->hexa[1])
+      printf(" %d", f->hexa[1]->Num);
+    if(f->prism[0])
+      printf(", in prism %d", f->prism[0]->Num);
+    if(f->prism[1])
+      printf(" %d", f->prism[1]->Num);
+    if(f->pyramid[0])
+      printf(", in pyramid %d", f->pyramid[0]->Num);
+    if(f->pyramid[1])
+      printf(" %d", f->pyramid[1]->Num);
+    printf("\n");
+  }
+  printf("Print QuadFacess STOP\n");
+  List_Delete(temp);
+}
+
+int compareQuadFace(const void *a, const void *b)
+{
+  QuadFace *q = (QuadFace *) a;
+  QuadFace *w = (QuadFace *) b;
+
+  if(q->V[0]->Num > w->V[0]->Num)
+    return 1;
+  if(q->V[0]->Num < w->V[0]->Num)
+    return -1;
+
+  if(q->V[1]->Num > w->V[1]->Num)
+    return 1;
+  if(q->V[1]->Num < w->V[1]->Num)
+    return -1;
+
+  if(q->V[2]->Num > w->V[2]->Num)
+    return 1;
+  if(q->V[2]->Num < w->V[2]->Num)
+    return -1;
+
+  if(q->V[3]->Num > w->V[3]->Num)
+    return 1;
+  if(q->V[3]->Num < w->V[3]->Num)
+    return -1;
+
+  return 0;
+}
diff --git a/Mesh/Face.h b/Mesh/Face.h
new file mode 100644
index 0000000000..1572a7e4c1
--- /dev/null
+++ b/Mesh/Face.h
@@ -0,0 +1,65 @@
+#ifndef _FACE_H_
+#define _FACE_H_
+
+// Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "List.h"
+#include "Tree.h"
+#include "Vertex.h"
+#include "Element.h"
+
+typedef struct {
+  Vertex *V[3];
+} Face;
+
+class QuadFace
+{
+ public:
+  Vertex *V[4];
+  Vertex *newv;
+  Quadrangle *quad;
+  Hexahedron *hexa[2];
+  Prism *prism[2];
+  Pyramid *pyramid[2];
+  QuadFace();
+  ~QuadFace(){;}
+};
+
+class QuadFacesContainer
+{
+ public :
+  Tree_T *AllFaces;
+  QuadFacesContainer();
+  ~QuadFacesContainer();
+  void AddQuadrangleTree(Tree_T *Quadrangles);
+  void AddHexahedronTree(Tree_T *Hexahedra);
+  void AddPrismTree(Tree_T *Prisms);    
+  void AddPyramidTree(Tree_T *Pyramids);
+  void AddQuadFaces(Quadrangle *q);
+  void AddQuadFaces(Hexahedron *h);
+  void AddQuadFaces(Prism *p);
+  void AddQuadFaces(Pyramid *p);
+  void Print();
+};
+
+int compareFace(const void *a, const void *b);
+int compareQuadFace(const void *a, const void *b);
+
+#endif
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 431832baef..da26bb7949 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.62 2004-05-25 04:10:05 geuzaine Exp $
+# $Id: Makefile,v 1.63 2004-05-25 23:16:26 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -67,6 +67,7 @@ SRC = 1D_Mesh.cpp \
       CrossData.cpp \
       Vertex.cpp \
       Edge.cpp \
+      Face.cpp \
       Element.cpp \
       Simplex.cpp 
 
@@ -100,228 +101,239 @@ depend:
 1D_Mesh.o: 1D_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
-  Mesh.h Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
-  Metric.h Matrix.h Utils.h ../Common/Context.h Interpolation.h
+  Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
+  ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h Utils.h \
+  ../Common/Context.h Interpolation.h
 2D_Mesh.o: 2D_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h Create.h 2D_Mesh.h \
-  ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h Create.h \
+  2D_Mesh.h ../Common/Context.h
 2D_SMesh.o: 2D_SMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h Mesh.h Vertex.h \
-  Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
+  Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
   Matrix.h ../Numeric/Numeric.h Interpolation.h
 2D_Elliptic.o: 2D_Elliptic.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h
 2D_BGMesh.o: 2D_BGMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h
 2D_Recombine.o: 2D_Recombine.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h Utils.h 2D_Mesh.h Create.h ../Common/Context.h
 2D_InitMesh.o: 2D_InitMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h
 2D_Bowyer.o: 2D_Bowyer.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h
 2D_Bricks.o: 2D_Bricks.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h
 2D_DivAndConq.o: 2D_DivAndConq.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h
 2D_Util.o: 2D_Util.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h ../Common/Context.h
 2D_Links.o: 2D_Links.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 2D_Mesh.h ../Common/Context.h
 2D_Tree.o: 2D_Tree.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
-  Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h 2D_Mesh.h
+  Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h 2D_Mesh.h
 2D_Cylindrical.o: 2D_Cylindrical.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h ../Common/Context.h
 2D_Parametric.o: 2D_Parametric.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h 2D_Mesh.h Create.h \
-  ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h \
+  2D_Mesh.h Create.h ../Common/Context.h
 2D_Mesh_Aniso.o: 2D_Mesh_Aniso.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h Create.h \
-  ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h \
+  Create.h ../Common/Context.h
 2D_Mesh_Triangle.o: 2D_Mesh_Triangle.cpp ../Common/Gmsh.h \
   ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h \
-  Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
+  Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
   Matrix.h ../Numeric/Numeric.h ../Common/Context.h
 3D_Mesh.o: 3D_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h 3D_Mesh.h Create.h ../Common/Context.h
 3D_SMesh.o: 3D_SMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
-  Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h Interpolation.h \
-  Create.h
+  Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h \
+  Interpolation.h Create.h
 3D_BGMesh.o: 3D_BGMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
-  Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h 2D_Mesh.h \
+  Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h 2D_Mesh.h \
   3D_Mesh.h ../Common/Views.h ../Common/ColorTable.h ../Numeric/Numeric.h \
   ../Common/Context.h
 3D_Extrude.o: 3D_Extrude.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Context.h Create.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Context.h \
+  Create.h
 3D_Extrude_Old.o: 3D_Extrude_Old.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Context.h Create.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Context.h \
+  Create.h
 3D_Coherence.o: 3D_Coherence.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
-  Mesh.h Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
-  Metric.h Matrix.h 3D_Mesh.h Create.h
+  Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
+  ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h 3D_Mesh.h Create.h
 3D_Divide.o: 3D_Divide.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h
 3D_Bricks.o: 3D_Bricks.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h
 MeshQuality.o: MeshQuality.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h
 Create.o: Create.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h ../Common/Context.h Create.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Utils.h \
+  ../Common/Context.h Create.h
 Generator.o: Generator.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h Create.h ../Common/Context.h ../Parser/OpenFile.h \
   ../Common/Views.h ../Common/ColorTable.h
 Print_Mesh.o: Print_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Create.h ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Create.h \
+  ../Common/Context.h
 Read_Mesh.o: Read_Mesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h ../Geo/CAD.h \
   ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Simplex.h \
-  ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h \
-  ../Mesh/Matrix.h 3D_Mesh.h Create.h ../Geo/MinMax.h ../Common/Context.h
+  ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
+  ../Mesh/Metric.h ../Mesh/Matrix.h 3D_Mesh.h Create.h ../Geo/MinMax.h \
+  ../Common/Context.h
 STL.o: STL.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
-  ../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Edge.h \
+  ../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
   ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h ../Geo/CAD.h \
   ../Geo/Geo.h Create.h Interpolation.h ../Common/Context.h
 SwapEdge.o: SwapEdge.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h SwapPatterns.h
 Utils.o: Utils.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h ../Common/Context.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h \
+  ../Common/Context.h
 Metric.o: Metric.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Interpolation.h
 Nurbs.o: Nurbs.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Nurbs.h Vertex.h Mesh.h Element.h \
-  Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h \
+  Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h \
   ../Geo/Geo.h ../Geo/GeoUtils.h Create.h ../Geo/CAD.h
 Interpolation.o: Interpolation.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
-  Nurbs.h Vertex.h Mesh.h Element.h Simplex.h Edge.h \
+  Nurbs.h Vertex.h Mesh.h Element.h Simplex.h Face.h Edge.h \
   ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h ../Geo/CAD.h Utils.h \
   Interpolation.h
 SecondOrder.o: SecondOrder.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Geo/Geo.h Mesh.h Vertex.h \
-  Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
+  Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
   Matrix.h Utils.h Interpolation.h ../Numeric/Numeric.h
 Smoothing.o: Smoothing.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Mesh.h \
-  Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
   Metric.h Matrix.h
 CrossData.o: CrossData.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
-  Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
+  Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
 Vertex.o: Vertex.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h Vertex.h \
-  Mesh.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h \
-  Matrix.h ../Common/Context.h
+  Mesh.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h STL.h \
+  Metric.h Matrix.h ../Common/Context.h
 Edge.o: Edge.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
-  ../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Edge.h \
+  ../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
+  ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
+Face.o: Face.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
+  ../Numeric/Numeric.h Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
   ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
 Element.o: Element.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h Mesh.h Vertex.h Element.h Simplex.h \
-  Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
+  Face.h Edge.h ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h
 Simplex.o: Simplex.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
-  Mesh.h Vertex.h Element.h Simplex.h Edge.h ../Geo/ExtrudeParams.h STL.h \
-  Metric.h Matrix.h ../Common/Context.h
+  Mesh.h Vertex.h Element.h Simplex.h Face.h Edge.h \
+  ../Geo/ExtrudeParams.h STL.h Metric.h Matrix.h ../Common/Context.h
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index ca05d083ea..56ca1379f9 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -488,10 +488,6 @@ void Freeze_Vertex(void *a, void *b);
 void deFreeze_Vertex(void *a, void *b);
 
 double Lc_XYZ(double X, double Y, double Z, Mesh * m);
-void Degre1();
-void Degre2(int dim);
-void Degre2_Curve(void *a, void *b);
-void Degre2_Surface(void *a, void *b);
 void ActionLiss(void *data, void *dummy);
 void ActionLissSurf(void *data, void *dummy);
 int  Recombine(Tree_T *TreeAllVert, Tree_T *TreeAllSimp, Tree_T *TreeAllQuad,
@@ -500,6 +496,11 @@ void ApplyLcFactor(Mesh *M);
 void ExportLcFieldOnVolume(Mesh * M, char *filename);
 void ExportLcFieldOnSurfaces(Mesh * M, char *filename);
 
+void Degre1();
+void Degre2(int dim);
+void Degre2_Curve(void *a, void *b);
+void Degre2_Surface(void *a, void *b);
+
 void Gamma_Maillage(Mesh * m, double *gamma, double *gammamax, double *gammamin);
 void Eta_Maillage(Mesh * m, double *gamma, double *gammamax, double *gammamin);
 void R_Maillage(Mesh * m, double *gamma, double *gammamax, double *gammamin);
diff --git a/Mesh/Print_Mesh.cpp b/Mesh/Print_Mesh.cpp
index 6e67ef23d1..31fd109a64 100644
--- a/Mesh/Print_Mesh.cpp
+++ b/Mesh/Print_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Print_Mesh.cpp,v 1.51 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Print_Mesh.cpp,v 1.52 2004-05-25 23:16:26 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -202,7 +202,7 @@ static void _msh_print_quadrangle(void *a, void *b)
   nbn = 4;
   if(q->VSUP) {
     type = QUADRANGLE_2;
-    nbs = 4;
+    nbs = 4 + 1;
   }
   else
     type = QUADRANGLE;
@@ -226,8 +226,11 @@ static void _msh_print_quadrangle(void *a, void *b)
   else {
     for(i = 0; i < nbn; i++)
       fprintf(MSHFILE, " %d", q->V[nbn - i - 1]->Num);
-    for(i = 0; i < nbs; i++)
-      fprintf(MSHFILE, " %d", q->VSUP[nbs - i - 1]->Num);
+    if(nbs){
+      for(i = 0; i < nbs - 1; i++)
+	fprintf(MSHFILE, " %d", q->VSUP[nbs - i - 2]->Num);
+      fprintf(MSHFILE, " %d", q->VSUP[nbs - 1]->Num);
+    }
   }
 
   fprintf(MSHFILE, "\n");
@@ -250,7 +253,7 @@ static void _msh_print_hexahedron(void *a, void *b)
   nbn = 8;
   if(h->VSUP) {
     type = HEXAHEDRON_2;
-    nbs = 12;
+    nbs = 12 + 6;
   }
   else
     type = HEXAHEDRON;
@@ -290,7 +293,7 @@ static void _msh_print_prism(void *a, void *b)
   nbn = 6;
   if(p->VSUP) {
     type = PRISM_2;
-    nbs = 9;
+    nbs = 9 + 3;
   }
   else {
     type = PRISM;
@@ -331,7 +334,7 @@ static void _msh_print_pyramid(void *a, void *b)
   nbn = 5;
   if(p->VSUP) {
     type = PYRAMID_2;
-    nbs = 8;
+    nbs = 8 + 1;
   }
   else {
     type = PYRAMID;
@@ -743,7 +746,7 @@ static void _unv_process_2D_elements(Mesh *m)
       for(int j = 0; j < List_Nbr(Elements); j++) {
         List_Read(Elements, j, &qx);
 	if(qx->VSUP)
-	  _unv_print_record(abs(qx->Num), QUAD, s->Num, 4, 4, &qx->V[0], qx->VSUP);
+	  _unv_print_record(abs(qx->Num), QUAD, s->Num, 4, 4+1, &qx->V[0], qx->VSUP);
 	else
 	  _unv_print_record(abs(qx->Num), QUAD2, s->Num, 4, 0, &qx->V[0], NULL);
       }
@@ -797,7 +800,7 @@ static void _unv_process_3D_elements(Mesh *m)
     for(int j = 0; j < List_Nbr(Elements); j++) {
       List_Read(Elements, j, &px);
       if(px->VSUP)
-	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 9, &px->V[0], px->VSUP);
+	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 9+3, &px->V[0], px->VSUP);
       else
 	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 0, &px->V[0], NULL);
     }
@@ -808,7 +811,7 @@ static void _unv_process_3D_elements(Mesh *m)
     for(int j = 0; j < List_Nbr(Elements); j++) {
       List_Read(Elements, j, &hx);
       if(hx->VSUP)
-	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 12, &hx->V[0], hx->VSUP);
+	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 12+6, &hx->V[0], hx->VSUP);
       else
 	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 0, &hx->V[0], NULL);
     }
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
index a53d0439d4..59fe5b1742 100644
--- a/Mesh/Read_Mesh.cpp
+++ b/Mesh/Read_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Read_Mesh.cpp,v 1.72 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Read_Mesh.cpp,v 1.73 2004-05-25 23:16:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -242,15 +242,15 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 	  case TRI1: Nbr_Nodes = 3; break;
 	  case TRI2: Nbr_Nodes = 6; break;
 	  case QUA1: Nbr_Nodes = 4; break;
-	  case QUA2: Nbr_Nodes = 8; break;
+	  case QUA2: Nbr_Nodes = 8 + 1; break;
 	  case TET1: Nbr_Nodes = 4; break;
 	  case TET2: Nbr_Nodes = 10; break;
 	  case HEX1: Nbr_Nodes = 8; break;
-	  case HEX2: Nbr_Nodes = 20; break;
+	  case HEX2: Nbr_Nodes = 20 + 6; break;
 	  case PRI1: Nbr_Nodes = 6; break;
-	  case PRI2: Nbr_Nodes = 15; break;
+	  case PRI2: Nbr_Nodes = 15 + 3; break;
 	  case PYR1: Nbr_Nodes = 5; break;
-	  case PYR2: Nbr_Nodes = 13; break;
+	  case PYR2: Nbr_Nodes = 13 + 1; break;
 	  }
 	}
 
@@ -330,8 +330,8 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           quad->iEnt = Elementary;
           quad->iPart = Add_MeshPartition(Partition, M);
 	  if(Type == QUA2){
-	    quad->VSUP = (Vertex **) Malloc(4 * sizeof(Vertex *));
-	    for(i = 0; i < 4; i++){
+	    quad->VSUP = (Vertex **) Malloc((4 + 1) * sizeof(Vertex *));
+	    for(i = 0; i < 4 + 1; i++){
 	      quad->VSUP[i] = vertsp[i+4];
 	      quad->VSUP[i]->Degree = 2;
 	    }
@@ -371,8 +371,8 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           hex->iEnt = Elementary;
           hex->iPart = Add_MeshPartition(Partition, M);
 	  if(Type == HEX2){
-	    hex->VSUP = (Vertex **) Malloc(12 * sizeof(Vertex *));
-	    for(i = 0; i < 12; i++){
+	    hex->VSUP = (Vertex **) Malloc((12 + 6) * sizeof(Vertex *));
+	    for(i = 0; i < 12 + 6; i++){
 	      hex->VSUP[i] = vertsp[i+8];
 	      hex->VSUP[i]->Degree = 2;
 	    }
@@ -392,8 +392,8 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           pri->iEnt = Elementary;
           pri->iPart = Add_MeshPartition(Partition, M);
 	  if(Type == PRI2){
-	    pri->VSUP = (Vertex **) Malloc(9 * sizeof(Vertex *));
-	    for(i = 0; i < 9; i++){
+	    pri->VSUP = (Vertex **) Malloc((9 + 3) * sizeof(Vertex *));
+	    for(i = 0; i < 9 + 3; i++){
 	      pri->VSUP[i] = vertsp[i+6];
 	      pri->VSUP[i]->Degree = 2;
 	    }
@@ -413,8 +413,8 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           pyr->iEnt = Elementary;
           pyr->iPart = Add_MeshPartition(Partition, M);
 	  if(Type == PYR2){
-	    pyr->VSUP = (Vertex **) Malloc(8 * sizeof(Vertex *));
-	    for(i = 0; i < 8; i++){
+	    pyr->VSUP = (Vertex **) Malloc((8 + 1) * sizeof(Vertex *));
+	    for(i = 0; i < 8 + 1; i++){
 	      pyr->VSUP[i] = vertsp[i+5];
 	      pyr->VSUP[i]->Degree = 2;
 	    }
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index 496d8eb231..e96239c2c9 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.23 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.24 2004-05-25 23:16:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -26,22 +26,24 @@
 #include "Interpolation.h"
 #include "Numeric.h"
 
-// FIXME: still to add middle face nodes for quads, hexas, prisms and
-// pyramids
-
 extern Mesh *THEM;
+
 extern int edges_tetra[6][2];
 extern int edges_quad[4][2];
 extern int edges_hexa[12][2];
 extern int edges_prism[9][2];
 extern int edges_pyramid[8][2];
 
+extern int quadfaces_hexa[6][4];
+extern int quadfaces_prism[3][4];
+
 static Surface *THES = NULL;
 static Curve *THEC = NULL;
 static EdgesContainer *edges = NULL;
+static QuadFacesContainer *faces = NULL;
 static List_T *VerticesToDelete = NULL;
 
-Vertex *oncurve(Vertex * v1, Vertex * v2)
+Vertex *onCurve(Vertex * v1, Vertex * v2)
 {
   Vertex v, *pv;
 
@@ -97,7 +99,7 @@ Vertex *oncurve(Vertex * v1, Vertex * v2)
   return pv;
 }
 
-Vertex *onsurface(Vertex * v1, Vertex * v2)
+Vertex *onSurface(Vertex * v1, Vertex * v2)
 {
   Vertex v, *pv;
   double U, V, U1, U2, V1, V2;
@@ -119,20 +121,43 @@ Vertex *onsurface(Vertex * v1, Vertex * v2)
   return pv;
 }
 
-void PutMiddlePoint(void *a, void *b)
+Vertex *onSurface(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
+{
+  Vertex v, *pv;
+  double U, V, U1, U2, U3, U4, V1, V2, V3, V4;
+
+  if(!THES)
+    return NULL;
+  if(THES->Typ == MSH_SURF_PLAN)
+    return NULL;
+
+  XYZtoUV(THES, v1->Pos.X, v1->Pos.Y, v1->Pos.Z, &U1, &V1, 1.0);
+  XYZtoUV(THES, v2->Pos.X, v2->Pos.Y, v2->Pos.Z, &U2, &V2, 1.0);
+  XYZtoUV(THES, v3->Pos.X, v3->Pos.Y, v3->Pos.Z, &U3, &V3, 1.0);
+  XYZtoUV(THES, v4->Pos.X, v4->Pos.Y, v4->Pos.Z, &U4, &V4, 1.0);
+
+  U = 0.25 * (U1 + U2 + U3 + U4);
+  V = 0.25 * (V1 + V2 + V3 + V4);
+  v = InterpolateSurface(THES, U, V, 0, 0);
+  pv = Create_Vertex(++THEM->MaxPointNum, v.Pos.X, v.Pos.Y, v.Pos.Z, v.lc, v.u);
+  pv->Degree = 2;
+
+  return pv;
+}
+
+void putMiddleEdgePoint(void *a, void *b)
 {
-  Edge *ed;
   Vertex *v;
-  int i, j, N;
+  int N;
 
-  ed = (Edge *) a;
+  Edge *ed = (Edge *) a;
 
   if(ed->newv){
     v = ed->newv;
   }
-  else if((v = oncurve(ed->V[0], ed->V[1]))){
+  else if((v = onCurve(ed->V[0], ed->V[1]))){
   }
-  else if((v = onsurface(ed->V[0], ed->V[1]))){
+  else if((v = onSurface(ed->V[0], ed->V[1]))){
   }
   else{
     v = Create_Vertex(++THEM->MaxPointNum,
@@ -147,7 +172,7 @@ void PutMiddlePoint(void *a, void *b)
   ed->newv = v;
   Tree_Insert(THEM->Vertices, &v);
       
-  for(i = 0; i < List_Nbr(ed->Simplexes); i++) {
+  for(int i = 0; i < List_Nbr(ed->Simplexes); i++) {
     Simplex *s;
     List_Read(ed->Simplexes, i, &s);
     if(s->V[3]) // tetrahedron
@@ -158,7 +183,7 @@ void PutMiddlePoint(void *a, void *b)
       N = 1;
     if(!s->VSUP)
       s->VSUP = (Vertex **) Malloc(N * sizeof(Vertex *));
-    for(j = 0; j < N; j++) {
+    for(int j = 0; j < N; j++) {
       if((!compareVertex(&s->V[edges_tetra[j][0]], &ed->V[0]) &&
 	  !compareVertex(&s->V[edges_tetra[j][1]], &ed->V[1])) ||
 	 (!compareVertex(&s->V[edges_tetra[j][0]], &ed->V[1]) &&
@@ -168,12 +193,12 @@ void PutMiddlePoint(void *a, void *b)
     }
   }
 
-  for(i = 0; i < List_Nbr(ed->Quadrangles); i++) {
+  for(int i = 0; i < List_Nbr(ed->Quadrangles); i++) {
     Quadrangle *q;
     List_Read(ed->Quadrangles, i, &q);
     if(!q->VSUP)
-      q->VSUP = (Vertex **) Malloc(4 * sizeof(Vertex *));
-    for(j = 0; j < 4; j++) {
+      q->VSUP = (Vertex **) Malloc(5 * sizeof(Vertex *));
+    for(int j = 0; j < 4; j++) {
       if((!compareVertex(&q->V[edges_quad[j][0]], &ed->V[0]) &&
           !compareVertex(&q->V[edges_quad[j][1]], &ed->V[1])) ||
          (!compareVertex(&q->V[edges_quad[j][0]], &ed->V[1]) &&
@@ -183,12 +208,12 @@ void PutMiddlePoint(void *a, void *b)
     }
   }
 
-  for(i = 0; i < List_Nbr(ed->Hexahedra); i++) {
+  for(int i = 0; i < List_Nbr(ed->Hexahedra); i++) {
     Hexahedron *h;
     List_Read(ed->Hexahedra, i, &h);
     if(!h->VSUP)
-      h->VSUP = (Vertex **) Malloc(12 * sizeof(Vertex *));
-    for(j = 0; j < 12; j++) {
+      h->VSUP = (Vertex **) Malloc(18 * sizeof(Vertex *));
+    for(int j = 0; j < 12; j++) {
       if((!compareVertex(&h->V[edges_hexa[j][0]], &ed->V[0]) &&
           !compareVertex(&h->V[edges_hexa[j][1]], &ed->V[1])) ||
          (!compareVertex(&h->V[edges_hexa[j][0]], &ed->V[1]) &&
@@ -198,12 +223,12 @@ void PutMiddlePoint(void *a, void *b)
     }
   }
 
-  for(i = 0; i < List_Nbr(ed->Prisms); i++) {
+  for(int i = 0; i < List_Nbr(ed->Prisms); i++) {
     Prism *p;
     List_Read(ed->Prisms, i, &p);
     if(!p->VSUP)
-      p->VSUP = (Vertex **) Malloc(9 * sizeof(Vertex *));
-    for(j = 0; j < 9; j++) {
+      p->VSUP = (Vertex **) Malloc(12 * sizeof(Vertex *));
+    for(int j = 0; j < 9; j++) {
       if((!compareVertex(&p->V[edges_prism[j][0]], &ed->V[0]) &&
           !compareVertex(&p->V[edges_prism[j][1]], &ed->V[1])) ||
          (!compareVertex(&p->V[edges_prism[j][0]], &ed->V[1]) &&
@@ -213,12 +238,12 @@ void PutMiddlePoint(void *a, void *b)
     }
   }
 
-  for(i = 0; i < List_Nbr(ed->Pyramids); i++) {
+  for(int i = 0; i < List_Nbr(ed->Pyramids); i++) {
     Pyramid *p;
     List_Read(ed->Pyramids, i, &p);
     if(!p->VSUP)
-      p->VSUP = (Vertex **) Malloc(8 * sizeof(Vertex *));
-    for(j = 0; j < 8; j++) {
+      p->VSUP = (Vertex **) Malloc(9 * sizeof(Vertex *));
+    for(int j = 0; j < 8; j++) {
       if((!compareVertex(&p->V[edges_pyramid[j][0]], &ed->V[0]) &&
           !compareVertex(&p->V[edges_pyramid[j][1]], &ed->V[1])) ||
          (!compareVertex(&p->V[edges_pyramid[j][0]], &ed->V[1]) &&
@@ -230,6 +255,81 @@ void PutMiddlePoint(void *a, void *b)
 
 }
 
+void putMiddleFacePoint(void *a, void *b)
+{
+  QuadFace *fac = (QuadFace*)a;
+  Vertex *v;
+
+  if(fac->newv){
+    v = fac->newv;
+  }
+  else if((v = onSurface(fac->V[0], fac->V[1], fac->V[2], fac->V[3]))){
+  }
+  else{
+    v = Create_Vertex(++THEM->MaxPointNum,
+                      0.25 * (fac->V[0]->Pos.X + fac->V[1]->Pos.X +
+			      fac->V[2]->Pos.X + fac->V[3]->Pos.X),
+                      0.25 * (fac->V[0]->Pos.Y + fac->V[1]->Pos.Y +
+			      fac->V[2]->Pos.Y + fac->V[3]->Pos.Y),
+                      0.25 * (fac->V[0]->Pos.Z + fac->V[1]->Pos.Z +
+			      fac->V[2]->Pos.Z + fac->V[3]->Pos.Z),
+                      0.25 * (fac->V[0]->lc + fac->V[1]->lc +
+			      fac->V[2]->lc + fac->V[3]->lc),
+                      0.25 * (fac->V[0]->u + fac->V[1]->u +
+			      fac->V[2]->u + fac->V[3]->u));
+    v->Degree = 2;
+  }
+
+  fac->newv = v;
+  Tree_Insert(THEM->Vertices, &v);
+
+  if(fac->quad && fac->quad->VSUP){
+    fac->quad->VSUP[4] = v;
+  }
+
+  QuadFace F;
+  
+  for(int i = 0; i < 2; i++) {
+    if(fac->hexa[i] && fac->hexa[i]->VSUP){
+      for(int j = 0; j < 6; j++) {
+	F.V[0] = fac->hexa[i]->V[quadfaces_hexa[j][0]];
+	F.V[1] = fac->hexa[i]->V[quadfaces_hexa[j][1]];
+	F.V[2] = fac->hexa[i]->V[quadfaces_hexa[j][2]];
+	F.V[3] = fac->hexa[i]->V[quadfaces_hexa[j][3]];
+	qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+	if(!compareQuadFace(fac, &F))
+	  fac->hexa[i]->VSUP[12+j] = v;
+      }
+    }
+  }
+
+  for(int i = 0; i < 2; i++) {
+    if(fac->prism[i] && fac->prism[i]->VSUP){
+      for(int j = 0; j < 3; j++) {
+	F.V[0] = fac->prism[i]->V[quadfaces_prism[j][0]];
+	F.V[1] = fac->prism[i]->V[quadfaces_prism[j][1]];
+	F.V[2] = fac->prism[i]->V[quadfaces_prism[j][2]];
+	F.V[3] = fac->prism[i]->V[quadfaces_prism[j][3]];
+	qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+	if(!compareQuadFace(fac, &F))
+	  fac->prism[i]->VSUP[9+j] = v;
+      }
+    }
+  }
+
+  for(int i = 0; i < 2; i++) {
+    if(fac->pyramid[i] && fac->pyramid[i]->VSUP){
+      F.V[0] = fac->pyramid[i]->V[0];
+      F.V[1] = fac->pyramid[i]->V[1];
+      F.V[2] = fac->pyramid[i]->V[2];
+      F.V[3] = fac->pyramid[i]->V[3];
+      qsort(F.V, 4, sizeof(Vertex *), compareVertex);
+      if(!compareQuadFace(fac, &F))
+	fac->pyramid[i]->VSUP[8] = v;
+    }
+  }
+}
+
 void ResetDegre2_Vertex(void *a, void *b)
 {
   Vertex *v = *(Vertex**)a;
@@ -299,11 +399,15 @@ void ResetDegre2_Volume(void *a, void *b)
 
 void Degre1()
 {
-  // (re-)initialize the global tree of edges
+  // (re-)initialize the global tree of edges/quadfaces
   if(edges)
     delete edges;
   edges = new EdgesContainer();
 
+  if(faces)
+    delete faces;
+  faces = new QuadFacesContainer();
+
   // reset VSUP in each element
   Tree_Action(THEM->Curves, ResetDegre2_Curve);
   Tree_Action(THEM->Surfaces, ResetDegre2_Surface);
@@ -325,21 +429,26 @@ void Degre2_Curve(void *a, void *b)
 {
   Curve *c = *(Curve**)a;
   if(c->Dirty) return;
+
   edges->AddSimplexTree(c->Simplexes);
   THEC = c;
   THES = NULL;
-  Tree_Action(edges->AllEdges, PutMiddlePoint);
+  Tree_Action(edges->AllEdges, putMiddleEdgePoint);
 }
 
 void Degre2_Surface(void *a, void *b)
 {
   Surface *s = *(Surface**)a;
   if(s->Dirty) return;
+
   edges->AddSimplexTree(s->Simplexes);
   edges->AddQuadrangleTree(s->Quadrangles);
   THEC = NULL;
   THES = s;
-  Tree_Action(edges->AllEdges, PutMiddlePoint);
+  Tree_Action(edges->AllEdges, putMiddleEdgePoint);
+
+  faces->AddQuadrangleTree(s->Quadrangles);
+  Tree_Action(faces->AllFaces, putMiddleFacePoint);
 }
 
 void Degre2_Volume(void *a, void *b)
@@ -353,7 +462,12 @@ void Degre2_Volume(void *a, void *b)
   edges->AddPyramidTree(v->Pyramids);
   THEC = NULL;
   THES = NULL;
-  Tree_Action(edges->AllEdges, PutMiddlePoint);
+  Tree_Action(edges->AllEdges, putMiddleEdgePoint);
+
+  faces->AddHexahedronTree(v->Hexahedra);
+  faces->AddPrismTree(v->Prisms);
+  faces->AddPyramidTree(v->Pyramids);
+  Tree_Action(faces->AllFaces, putMiddleFacePoint);
 }
 
 void Degre2(int dim)
diff --git a/Mesh/Simplex.cpp b/Mesh/Simplex.cpp
index f31ce24382..a34f667420 100644
--- a/Mesh/Simplex.cpp
+++ b/Mesh/Simplex.cpp
@@ -1,4 +1,4 @@
-// $Id: Simplex.cpp,v 1.32 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Simplex.cpp,v 1.33 2004-05-25 23:16:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -782,30 +782,3 @@ bool Simplex::SwapFace(int iFac, List_T * newsimp, List_T * delsimp)
   return true;
 }
 
-int compareFace(const void *a, const void *b)
-{
-  Face *q = (Face *) a;
-  Face *w = (Face *) b;
-
-  if(!q->V[0] || !w->V[0])
-    Msg(FATAL, "Bad face (are you trying to generate hexahedra with a Delaunay?!)");
-
-  if(q->V[0]->Num > w->V[0]->Num)
-    return (1);
-  if(q->V[0]->Num < w->V[0]->Num)
-    return (-1);
-
-  if(q->V[1]->Num > w->V[1]->Num)
-    return (1);
-  if(q->V[1]->Num < w->V[1]->Num)
-    return (-1);
-
-  if(FACE_DIMENSION == 1 || !q->V[2] || !w->V[2])
-    return 0;
-
-  if(q->V[2]->Num > w->V[2]->Num)
-    return (1);
-  if(q->V[2]->Num < w->V[2]->Num)
-    return (-1);
-  return (0);
-}
diff --git a/Mesh/Simplex.h b/Mesh/Simplex.h
index 451c039e06..167e4ba814 100644
--- a/Mesh/Simplex.h
+++ b/Mesh/Simplex.h
@@ -23,10 +23,7 @@
 #include "List.h"
 #include "Vertex.h"
 #include "Element.h"
-
-typedef struct {
-  Vertex *V[3];
-}Face;
+#include "Face.h"
 
 class Simplex : public Element {
  public:
diff --git a/Mesh/SwapEdge.cpp b/Mesh/SwapEdge.cpp
index 263b043815..47b96bb612 100644
--- a/Mesh/SwapEdge.cpp
+++ b/Mesh/SwapEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: SwapEdge.cpp,v 1.15 2004-02-07 01:40:22 geuzaine Exp $
+// $Id: SwapEdge.cpp,v 1.16 2004-05-25 23:16:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -65,7 +65,7 @@ int TrouveCoquille(Simplex * s, Vertex * Ed[2],
   N = 0;
   Contour[N++] = E.V[0];
   Contour[N++] = E.V[1];
-  List_Suppress(Edges, &E, compareedge);
+  List_Suppress(Edges, &E, compareEdge);
   K = 0;
 
   while(List_Nbr(Edges)) {
@@ -76,13 +76,13 @@ int TrouveCoquille(Simplex * s, Vertex * Ed[2],
       if(!compareVertex(&Contour[N - 1], &E.V[0]) &&
          compareVertex(&Contour[N - 2], &E.V[1])) {
         Contour[N++] = E.V[1];
-        List_Suppress(Edges, &E, compareedge);
+        List_Suppress(Edges, &E, compareEdge);
         break;
       }
       if(!compareVertex(&Contour[N - 1], &E.V[1]) &&
          compareVertex(&Contour[N - 2], &E.V[0])) {
         Contour[N++] = E.V[0];
-        List_Suppress(Edges, &E, compareedge);
+        List_Suppress(Edges, &E, compareEdge);
         break;
       }
     }
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 017b1a8688..41d53a4d29 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -191,7 +191,7 @@
 
 #line 1 "Gmsh.y"
 
-// $Id: Gmsh.tab.cpp,v 1.191 2004-05-25 04:10:05 geuzaine Exp $
+// $Id: Gmsh.tab.cpp,v 1.192 2004-05-25 23:16:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -757,21 +757,21 @@ static const short yyrline[] = { 0,
   2678,  2689,  2694,  2705,  2710,  2721,  2726,  2739,  2757,  2775,
   2793,  2798,  2816,  2821,  2839,  2844,  2865,  2882,  2899,  2916,
   2921,  2938,  2944,  2961,  2967,  2986,  2990,  2995,  3022,  3046,
-  3054,  3073,  3091,  3109,  3136,  3162,  3188,  3202,  3221,  3226,
-  3235,  3237,  3238,  3239,  3240,  3243,  3245,  3246,  3247,  3248,
-  3249,  3250,  3251,  3252,  3259,  3260,  3261,  3262,  3263,  3264,
-  3265,  3266,  3267,  3268,  3269,  3270,  3271,  3272,  3273,  3274,
-  3275,  3276,  3277,  3278,  3279,  3280,  3281,  3282,  3283,  3284,
-  3285,  3286,  3287,  3288,  3289,  3290,  3292,  3293,  3294,  3295,
-  3296,  3297,  3298,  3299,  3300,  3301,  3302,  3303,  3304,  3305,
-  3306,  3307,  3308,  3309,  3310,  3311,  3312,  3317,  3322,  3323,
-  3324,  3328,  3340,  3359,  3372,  3384,  3406,  3423,  3440,  3457,
-  3476,  3481,  3485,  3489,  3493,  3499,  3504,  3508,  3512,  3518,
-  3522,  3527,  3531,  3536,  3540,  3544,  3550,  3556,  3563,  3569,
-  3573,  3577,  3588,  3595,  3606,  3626,  3636,  3646,  3658,  3674,
-  3692,  3715,  3742,  3748,  3752,  3756,  3768,  3773,  3785,  3791,
-  3811,  3816,  3829,  3835,  3841,  3846,  3854,  3868,  3872,  3891,
-  3907
+  3054,  3073,  3091,  3109,  3136,  3162,  3188,  3202,  3220,  3225,
+  3234,  3236,  3237,  3238,  3239,  3242,  3244,  3245,  3246,  3247,
+  3248,  3249,  3250,  3251,  3258,  3259,  3260,  3261,  3262,  3263,
+  3264,  3265,  3266,  3267,  3268,  3269,  3270,  3271,  3272,  3273,
+  3274,  3275,  3276,  3277,  3278,  3279,  3280,  3281,  3282,  3283,
+  3284,  3285,  3286,  3287,  3288,  3289,  3291,  3292,  3293,  3294,
+  3295,  3296,  3297,  3298,  3299,  3300,  3301,  3302,  3303,  3304,
+  3305,  3306,  3307,  3308,  3309,  3310,  3311,  3316,  3321,  3322,
+  3323,  3327,  3339,  3358,  3371,  3383,  3405,  3422,  3439,  3456,
+  3475,  3480,  3484,  3488,  3492,  3498,  3503,  3507,  3511,  3517,
+  3521,  3526,  3530,  3535,  3539,  3543,  3549,  3555,  3562,  3568,
+  3572,  3576,  3587,  3594,  3605,  3625,  3635,  3645,  3657,  3673,
+  3691,  3714,  3741,  3747,  3751,  3755,  3767,  3772,  3784,  3790,
+  3810,  3815,  3828,  3834,  3840,  3845,  3853,  3867,  3871,  3890,
+  3906
 };
 #endif
 
@@ -7844,78 +7844,77 @@ case 338:
         Surface *s = FindSurface(j, THEM);
 	if(s){
 	  s->Recombine = 1;
-	  s->RecombineAngle = 30.;
         }
       }
       List_Delete(yyvsp[-1].l);
     ;
     break;}
 case 339:
-#line 3223 "Gmsh.y"
+#line 3222 "Gmsh.y"
 { 
       ReplaceAllDuplicates(THEM);
     ;
     break;}
 case 340:
-#line 3227 "Gmsh.y"
+#line 3226 "Gmsh.y"
 { 
       IntersectAllSegmentsTogether();
     ;
     break;}
 case 341:
-#line 3236 "Gmsh.y"
+#line 3235 "Gmsh.y"
 {yyval.i = 1;;
     break;}
 case 342:
-#line 3237 "Gmsh.y"
+#line 3236 "Gmsh.y"
 {yyval.i = 0;;
     break;}
 case 343:
-#line 3238 "Gmsh.y"
+#line 3237 "Gmsh.y"
 {yyval.i = -1;;
     break;}
 case 344:
-#line 3239 "Gmsh.y"
+#line 3238 "Gmsh.y"
 {yyval.i = -1;;
     break;}
 case 345:
-#line 3240 "Gmsh.y"
+#line 3239 "Gmsh.y"
 {yyval.i = -1;;
     break;}
 case 346:
-#line 3244 "Gmsh.y"
+#line 3243 "Gmsh.y"
 { yyval.d = yyvsp[0].d;           ;
     break;}
 case 347:
-#line 3245 "Gmsh.y"
+#line 3244 "Gmsh.y"
 { yyval.d = yyvsp[-1].d;           ;
     break;}
 case 348:
-#line 3246 "Gmsh.y"
+#line 3245 "Gmsh.y"
 { yyval.d = -yyvsp[0].d;          ;
     break;}
 case 349:
-#line 3247 "Gmsh.y"
+#line 3246 "Gmsh.y"
 { yyval.d = yyvsp[0].d;           ;
     break;}
 case 350:
-#line 3248 "Gmsh.y"
+#line 3247 "Gmsh.y"
 { yyval.d = !yyvsp[0].d;          ;
     break;}
 case 351:
-#line 3249 "Gmsh.y"
+#line 3248 "Gmsh.y"
 { yyval.d = yyvsp[-2].d - yyvsp[0].d;      ;
     break;}
 case 352:
-#line 3250 "Gmsh.y"
+#line 3249 "Gmsh.y"
 { yyval.d = yyvsp[-2].d + yyvsp[0].d;      ;
     break;}
 case 353:
-#line 3251 "Gmsh.y"
+#line 3250 "Gmsh.y"
 { yyval.d = yyvsp[-2].d * yyvsp[0].d;      ;
     break;}
 case 354:
-#line 3253 "Gmsh.y"
+#line 3252 "Gmsh.y"
 { 
       if(!yyvsp[0].d)
 	yymsg(GERROR, "Division by zero in '%g / %g'", yyvsp[-2].d, yyvsp[0].d);
@@ -7924,235 +7923,235 @@ case 354:
     ;
     break;}
 case 355:
-#line 3259 "Gmsh.y"
+#line 3258 "Gmsh.y"
 { yyval.d = (int)yyvsp[-2].d % (int)yyvsp[0].d;  ;
     break;}
 case 356:
-#line 3260 "Gmsh.y"
+#line 3259 "Gmsh.y"
 { yyval.d = pow(yyvsp[-2].d, yyvsp[0].d);  ;
     break;}
 case 357:
-#line 3261 "Gmsh.y"
+#line 3260 "Gmsh.y"
 { yyval.d = yyvsp[-2].d < yyvsp[0].d;      ;
     break;}
 case 358:
-#line 3262 "Gmsh.y"
+#line 3261 "Gmsh.y"
 { yyval.d = yyvsp[-2].d > yyvsp[0].d;      ;
     break;}
 case 359:
-#line 3263 "Gmsh.y"
+#line 3262 "Gmsh.y"
 { yyval.d = yyvsp[-2].d <= yyvsp[0].d;     ;
     break;}
 case 360:
-#line 3264 "Gmsh.y"
+#line 3263 "Gmsh.y"
 { yyval.d = yyvsp[-2].d >= yyvsp[0].d;     ;
     break;}
 case 361:
-#line 3265 "Gmsh.y"
+#line 3264 "Gmsh.y"
 { yyval.d = yyvsp[-2].d == yyvsp[0].d;     ;
     break;}
 case 362:
-#line 3266 "Gmsh.y"
+#line 3265 "Gmsh.y"
 { yyval.d = yyvsp[-2].d != yyvsp[0].d;     ;
     break;}
 case 363:
-#line 3267 "Gmsh.y"
+#line 3266 "Gmsh.y"
 { yyval.d = yyvsp[-2].d && yyvsp[0].d;     ;
     break;}
 case 364:
-#line 3268 "Gmsh.y"
+#line 3267 "Gmsh.y"
 { yyval.d = yyvsp[-2].d || yyvsp[0].d;     ;
     break;}
 case 365:
-#line 3269 "Gmsh.y"
+#line 3268 "Gmsh.y"
 { yyval.d = yyvsp[-4].d? yyvsp[-2].d : yyvsp[0].d;  ;
     break;}
 case 366:
-#line 3270 "Gmsh.y"
+#line 3269 "Gmsh.y"
 { yyval.d = exp(yyvsp[-1].d);      ;
     break;}
 case 367:
-#line 3271 "Gmsh.y"
+#line 3270 "Gmsh.y"
 { yyval.d = log(yyvsp[-1].d);      ;
     break;}
 case 368:
-#line 3272 "Gmsh.y"
+#line 3271 "Gmsh.y"
 { yyval.d = log10(yyvsp[-1].d);    ;
     break;}
 case 369:
-#line 3273 "Gmsh.y"
+#line 3272 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-1].d);     ;
     break;}
 case 370:
-#line 3274 "Gmsh.y"
+#line 3273 "Gmsh.y"
 { yyval.d = sin(yyvsp[-1].d);      ;
     break;}
 case 371:
-#line 3275 "Gmsh.y"
+#line 3274 "Gmsh.y"
 { yyval.d = asin(yyvsp[-1].d);     ;
     break;}
 case 372:
-#line 3276 "Gmsh.y"
+#line 3275 "Gmsh.y"
 { yyval.d = cos(yyvsp[-1].d);      ;
     break;}
 case 373:
-#line 3277 "Gmsh.y"
+#line 3276 "Gmsh.y"
 { yyval.d = acos(yyvsp[-1].d);     ;
     break;}
 case 374:
-#line 3278 "Gmsh.y"
+#line 3277 "Gmsh.y"
 { yyval.d = tan(yyvsp[-1].d);      ;
     break;}
 case 375:
-#line 3279 "Gmsh.y"
+#line 3278 "Gmsh.y"
 { yyval.d = atan(yyvsp[-1].d);     ;
     break;}
 case 376:
-#line 3280 "Gmsh.y"
+#line 3279 "Gmsh.y"
 { yyval.d = atan2(yyvsp[-3].d, yyvsp[-1].d);;
     break;}
 case 377:
-#line 3281 "Gmsh.y"
+#line 3280 "Gmsh.y"
 { yyval.d = sinh(yyvsp[-1].d);     ;
     break;}
 case 378:
-#line 3282 "Gmsh.y"
+#line 3281 "Gmsh.y"
 { yyval.d = cosh(yyvsp[-1].d);     ;
     break;}
 case 379:
-#line 3283 "Gmsh.y"
+#line 3282 "Gmsh.y"
 { yyval.d = tanh(yyvsp[-1].d);     ;
     break;}
 case 380:
-#line 3284 "Gmsh.y"
+#line 3283 "Gmsh.y"
 { yyval.d = fabs(yyvsp[-1].d);     ;
     break;}
 case 381:
-#line 3285 "Gmsh.y"
+#line 3284 "Gmsh.y"
 { yyval.d = floor(yyvsp[-1].d);    ;
     break;}
 case 382:
-#line 3286 "Gmsh.y"
+#line 3285 "Gmsh.y"
 { yyval.d = ceil(yyvsp[-1].d);     ;
     break;}
 case 383:
-#line 3287 "Gmsh.y"
+#line 3286 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 384:
-#line 3288 "Gmsh.y"
+#line 3287 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 385:
-#line 3289 "Gmsh.y"
+#line 3288 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-3].d*yyvsp[-3].d+yyvsp[-1].d*yyvsp[-1].d); ;
     break;}
 case 386:
-#line 3290 "Gmsh.y"
+#line 3289 "Gmsh.y"
 { yyval.d = yyvsp[-1].d*(double)rand()/(double)RAND_MAX; ;
     break;}
 case 387:
-#line 3292 "Gmsh.y"
+#line 3291 "Gmsh.y"
 { yyval.d = exp(yyvsp[-1].d);      ;
     break;}
 case 388:
-#line 3293 "Gmsh.y"
+#line 3292 "Gmsh.y"
 { yyval.d = log(yyvsp[-1].d);      ;
     break;}
 case 389:
-#line 3294 "Gmsh.y"
+#line 3293 "Gmsh.y"
 { yyval.d = log10(yyvsp[-1].d);    ;
     break;}
 case 390:
-#line 3295 "Gmsh.y"
+#line 3294 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-1].d);     ;
     break;}
 case 391:
-#line 3296 "Gmsh.y"
+#line 3295 "Gmsh.y"
 { yyval.d = sin(yyvsp[-1].d);      ;
     break;}
 case 392:
-#line 3297 "Gmsh.y"
+#line 3296 "Gmsh.y"
 { yyval.d = asin(yyvsp[-1].d);     ;
     break;}
 case 393:
-#line 3298 "Gmsh.y"
+#line 3297 "Gmsh.y"
 { yyval.d = cos(yyvsp[-1].d);      ;
     break;}
 case 394:
-#line 3299 "Gmsh.y"
+#line 3298 "Gmsh.y"
 { yyval.d = acos(yyvsp[-1].d);     ;
     break;}
 case 395:
-#line 3300 "Gmsh.y"
+#line 3299 "Gmsh.y"
 { yyval.d = tan(yyvsp[-1].d);      ;
     break;}
 case 396:
-#line 3301 "Gmsh.y"
+#line 3300 "Gmsh.y"
 { yyval.d = atan(yyvsp[-1].d);     ;
     break;}
 case 397:
-#line 3302 "Gmsh.y"
+#line 3301 "Gmsh.y"
 { yyval.d = atan2(yyvsp[-3].d, yyvsp[-1].d);;
     break;}
 case 398:
-#line 3303 "Gmsh.y"
+#line 3302 "Gmsh.y"
 { yyval.d = sinh(yyvsp[-1].d);     ;
     break;}
 case 399:
-#line 3304 "Gmsh.y"
+#line 3303 "Gmsh.y"
 { yyval.d = cosh(yyvsp[-1].d);     ;
     break;}
 case 400:
-#line 3305 "Gmsh.y"
+#line 3304 "Gmsh.y"
 { yyval.d = tanh(yyvsp[-1].d);     ;
     break;}
 case 401:
-#line 3306 "Gmsh.y"
+#line 3305 "Gmsh.y"
 { yyval.d = fabs(yyvsp[-1].d);     ;
     break;}
 case 402:
-#line 3307 "Gmsh.y"
+#line 3306 "Gmsh.y"
 { yyval.d = floor(yyvsp[-1].d);    ;
     break;}
 case 403:
-#line 3308 "Gmsh.y"
+#line 3307 "Gmsh.y"
 { yyval.d = ceil(yyvsp[-1].d);     ;
     break;}
 case 404:
-#line 3309 "Gmsh.y"
+#line 3308 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 405:
-#line 3310 "Gmsh.y"
+#line 3309 "Gmsh.y"
 { yyval.d = fmod(yyvsp[-3].d, yyvsp[-1].d); ;
     break;}
 case 406:
-#line 3311 "Gmsh.y"
+#line 3310 "Gmsh.y"
 { yyval.d = sqrt(yyvsp[-3].d*yyvsp[-3].d+yyvsp[-1].d*yyvsp[-1].d); ;
     break;}
 case 407:
-#line 3312 "Gmsh.y"
+#line 3311 "Gmsh.y"
 { yyval.d = yyvsp[-1].d*(double)rand()/(double)RAND_MAX; ;
     break;}
 case 408:
-#line 3321 "Gmsh.y"
+#line 3320 "Gmsh.y"
 { yyval.d = yyvsp[0].d; ;
     break;}
 case 409:
-#line 3322 "Gmsh.y"
+#line 3321 "Gmsh.y"
 { yyval.d = 3.141592653589793; ;
     break;}
 case 410:
-#line 3323 "Gmsh.y"
+#line 3322 "Gmsh.y"
 { yyval.d = ParUtil::Instance()->rank(); ;
     break;}
 case 411:
-#line 3324 "Gmsh.y"
+#line 3323 "Gmsh.y"
 { yyval.d = ParUtil::Instance()->size(); ;
     break;}
 case 412:
-#line 3329 "Gmsh.y"
+#line 3328 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[0].c;
@@ -8166,7 +8165,7 @@ case 412:
     ;
     break;}
 case 413:
-#line 3341 "Gmsh.y"
+#line 3340 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-3].c;
@@ -8187,7 +8186,7 @@ case 413:
     ;
     break;}
 case 414:
-#line 3360 "Gmsh.y"
+#line 3359 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-2].c;
@@ -8202,7 +8201,7 @@ case 414:
     ;
     break;}
 case 415:
-#line 3373 "Gmsh.y"
+#line 3372 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-1].c;
@@ -8216,7 +8215,7 @@ case 415:
     ;
     break;}
 case 416:
-#line 3385 "Gmsh.y"
+#line 3384 "Gmsh.y"
 {
       Symbol TheSymbol;
       TheSymbol.Name = yyvsp[-4].c;
@@ -8237,7 +8236,7 @@ case 416:
     ;
     break;}
 case 417:
-#line 3407 "Gmsh.y"
+#line 3406 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -8256,7 +8255,7 @@ case 417:
     ;
     break;}
 case 418:
-#line 3424 "Gmsh.y"
+#line 3423 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -8275,7 +8274,7 @@ case 418:
     ;
     break;}
 case 419:
-#line 3441 "Gmsh.y"
+#line 3440 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -8294,7 +8293,7 @@ case 419:
     ;
     break;}
 case 420:
-#line 3458 "Gmsh.y"
+#line 3457 "Gmsh.y"
 {
       double (*pNumOpt)(int num, int action, double value);
       StringXNumber *pNumCat;
@@ -8313,130 +8312,130 @@ case 420:
     ;
     break;}
 case 421:
-#line 3478 "Gmsh.y"
+#line 3477 "Gmsh.y"
 {
       memcpy(yyval.v, yyvsp[0].v, 5*sizeof(double));
     ;
     break;}
 case 422:
-#line 3482 "Gmsh.y"
+#line 3481 "Gmsh.y"
 {
       for(int i = 0; i < 5; i++) yyval.v[i] = -yyvsp[0].v[i];
     ;
     break;}
 case 423:
-#line 3486 "Gmsh.y"
+#line 3485 "Gmsh.y"
 { 
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[0].v[i];
     ;
     break;}
 case 424:
-#line 3490 "Gmsh.y"
+#line 3489 "Gmsh.y"
 { 
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[-2].v[i] - yyvsp[0].v[i];
     ;
     break;}
 case 425:
-#line 3494 "Gmsh.y"
+#line 3493 "Gmsh.y"
 {
       for(int i = 0; i < 5; i++) yyval.v[i] = yyvsp[-2].v[i] + yyvsp[0].v[i];
     ;
     break;}
 case 426:
-#line 3501 "Gmsh.y"
+#line 3500 "Gmsh.y"
 { 
       yyval.v[0] = yyvsp[-9].d;  yyval.v[1] = yyvsp[-7].d;  yyval.v[2] = yyvsp[-5].d;  yyval.v[3] = yyvsp[-3].d; yyval.v[4] = yyvsp[-1].d;
     ;
     break;}
 case 427:
-#line 3505 "Gmsh.y"
+#line 3504 "Gmsh.y"
 { 
       yyval.v[0] = yyvsp[-7].d;  yyval.v[1] = yyvsp[-5].d;  yyval.v[2] = yyvsp[-3].d;  yyval.v[3] = yyvsp[-1].d; yyval.v[4] = 1.0;
     ;
     break;}
 case 428:
-#line 3509 "Gmsh.y"
+#line 3508 "Gmsh.y"
 {
       yyval.v[0] = yyvsp[-5].d;  yyval.v[1] = yyvsp[-3].d;  yyval.v[2] = yyvsp[-1].d;  yyval.v[3] = 0.0; yyval.v[4] = 1.0;
     ;
     break;}
 case 429:
-#line 3513 "Gmsh.y"
+#line 3512 "Gmsh.y"
 {
       yyval.v[0] = yyvsp[-5].d;  yyval.v[1] = yyvsp[-3].d;  yyval.v[2] = yyvsp[-1].d;  yyval.v[3] = 0.0; yyval.v[4] = 1.0;
     ;
     break;}
 case 430:
-#line 3520 "Gmsh.y"
+#line 3519 "Gmsh.y"
 {
     ;
     break;}
 case 431:
-#line 3523 "Gmsh.y"
+#line 3522 "Gmsh.y"
 {
     ;
     break;}
 case 432:
-#line 3529 "Gmsh.y"
+#line 3528 "Gmsh.y"
 {
     ;
     break;}
 case 433:
-#line 3532 "Gmsh.y"
+#line 3531 "Gmsh.y"
 {
     ;
     break;}
 case 434:
-#line 3538 "Gmsh.y"
+#line 3537 "Gmsh.y"
 {
     ;
     break;}
 case 435:
-#line 3541 "Gmsh.y"
+#line 3540 "Gmsh.y"
 {
        yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 436:
-#line 3545 "Gmsh.y"
+#line 3544 "Gmsh.y"
 {
        yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 437:
-#line 3552 "Gmsh.y"
+#line 3551 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(List_T*));
       List_Add(yyval.l, &(yyvsp[0].l));
     ;
     break;}
 case 438:
-#line 3557 "Gmsh.y"
+#line 3556 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].l));
     ;
     break;}
 case 439:
-#line 3565 "Gmsh.y"
+#line 3564 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 440:
-#line 3570 "Gmsh.y"
+#line 3569 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
     ;
     break;}
 case 441:
-#line 3574 "Gmsh.y"
+#line 3573 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 442:
-#line 3578 "Gmsh.y"
+#line 3577 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
       double *pd;
@@ -8447,7 +8446,7 @@ case 442:
     ;
     break;}
 case 443:
-#line 3590 "Gmsh.y"
+#line 3589 "Gmsh.y"
 { 
       yyval.l = List_Create(2, 1, sizeof(double)); 
       for(double d = yyvsp[-2].d; (yyvsp[-2].d < yyvsp[0].d) ? (d <= yyvsp[0].d) : (d >= yyvsp[0].d); (yyvsp[-2].d < yyvsp[0].d) ? (d += 1.) : (d -= 1.)) 
@@ -8455,7 +8454,7 @@ case 443:
     ;
     break;}
 case 444:
-#line 3596 "Gmsh.y"
+#line 3595 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double)); 
       if(!yyvsp[0].d || (yyvsp[-4].d < yyvsp[-2].d && yyvsp[0].d < 0) || (yyvsp[-4].d > yyvsp[-2].d && yyvsp[0].d > 0)){
@@ -8468,7 +8467,7 @@ case 444:
    ;
     break;}
 case 445:
-#line 3607 "Gmsh.y"
+#line 3606 "Gmsh.y"
 {
       // Returns the coordinates of a point and fills a list with it.
       // This allows to ensure e.g. that relative point positions are
@@ -8490,7 +8489,7 @@ case 445:
     ;
     break;}
 case 446:
-#line 3627 "Gmsh.y"
+#line 3626 "Gmsh.y"
 {
       yyval.l = List_Create(List_Nbr(yyvsp[0].l), 1, sizeof(double));
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
@@ -8502,7 +8501,7 @@ case 446:
     ;
     break;}
 case 447:
-#line 3637 "Gmsh.y"
+#line 3636 "Gmsh.y"
 {
       yyval.l = List_Create(List_Nbr(yyvsp[0].l), 1, sizeof(double));
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
@@ -8514,7 +8513,7 @@ case 447:
     ;
     break;}
 case 448:
-#line 3647 "Gmsh.y"
+#line 3646 "Gmsh.y"
 {
       // FIXME: The syntax for this is ugly: we get double semi-colons
       // at the end of the line
@@ -8528,7 +8527,7 @@ case 448:
     ;
     break;}
 case 449:
-#line 3659 "Gmsh.y"
+#line 3658 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -8546,7 +8545,7 @@ case 449:
     ;
     break;}
 case 450:
-#line 3675 "Gmsh.y"
+#line 3674 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -8566,7 +8565,7 @@ case 450:
     ;
     break;}
 case 451:
-#line 3693 "Gmsh.y"
+#line 3692 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -8591,7 +8590,7 @@ case 451:
     ;
     break;}
 case 452:
-#line 3716 "Gmsh.y"
+#line 3715 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       Symbol TheSymbol;
@@ -8618,26 +8617,26 @@ case 452:
     ;
     break;}
 case 453:
-#line 3744 "Gmsh.y"
+#line 3743 "Gmsh.y"
 {
       yyval.l = List_Create(2, 1, sizeof(double));
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 454:
-#line 3749 "Gmsh.y"
+#line 3748 "Gmsh.y"
 {
       yyval.l = yyvsp[0].l;
     ;
     break;}
 case 455:
-#line 3753 "Gmsh.y"
+#line 3752 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].d));
     ;
     break;}
 case 456:
-#line 3757 "Gmsh.y"
+#line 3756 "Gmsh.y"
 {
       for(int i = 0; i < List_Nbr(yyvsp[0].l); i++){
 	double d;
@@ -8648,19 +8647,19 @@ case 456:
     ;
     break;}
 case 457:
-#line 3770 "Gmsh.y"
+#line 3769 "Gmsh.y"
 {
       yyval.u = PACK_COLOR((int)yyvsp[-7].d, (int)yyvsp[-5].d, (int)yyvsp[-3].d, (int)yyvsp[-1].d);
     ;
     break;}
 case 458:
-#line 3774 "Gmsh.y"
+#line 3773 "Gmsh.y"
 {
       yyval.u = PACK_COLOR((int)yyvsp[-5].d, (int)yyvsp[-3].d, (int)yyvsp[-1].d, 255);
     ;
     break;}
 case 459:
-#line 3786 "Gmsh.y"
+#line 3785 "Gmsh.y"
 {
       int flag;
       yyval.u = Get_ColorForString(ColorString, -1, yyvsp[0].c, &flag);
@@ -8668,7 +8667,7 @@ case 459:
     ;
     break;}
 case 460:
-#line 3792 "Gmsh.y"
+#line 3791 "Gmsh.y"
 {
       unsigned int (*pColOpt)(int num, int action, unsigned int value);
       StringXColor *pColCat;
@@ -8688,13 +8687,13 @@ case 460:
     ;
     break;}
 case 461:
-#line 3813 "Gmsh.y"
+#line 3812 "Gmsh.y"
 {
       yyval.l = yyvsp[-1].l;
     ;
     break;}
 case 462:
-#line 3817 "Gmsh.y"
+#line 3816 "Gmsh.y"
 {
       yyval.l = List_Create(256, 10, sizeof(unsigned int));
       GmshColorTable *ct = Get_ColorTable((int)yyvsp[-3].d);
@@ -8707,26 +8706,26 @@ case 462:
     ;
     break;}
 case 463:
-#line 3831 "Gmsh.y"
+#line 3830 "Gmsh.y"
 {
       yyval.l = List_Create(256, 10, sizeof(unsigned int));
       List_Add(yyval.l, &(yyvsp[0].u));
     ;
     break;}
 case 464:
-#line 3836 "Gmsh.y"
+#line 3835 "Gmsh.y"
 {
       List_Add(yyval.l, &(yyvsp[0].u));
     ;
     break;}
 case 465:
-#line 3843 "Gmsh.y"
+#line 3842 "Gmsh.y"
 {
       yyval.c = yyvsp[0].c;
     ;
     break;}
 case 466:
-#line 3847 "Gmsh.y"
+#line 3846 "Gmsh.y"
 {
       yyval.c = (char *)Malloc((strlen(yyvsp[-3].c)+strlen(yyvsp[-1].c)+1)*sizeof(char));
       strcpy(yyval.c, yyvsp[-3].c);  
@@ -8736,7 +8735,7 @@ case 466:
     ;
     break;}
 case 467:
-#line 3855 "Gmsh.y"
+#line 3854 "Gmsh.y"
 {
       yyval.c = (char *)Malloc((strlen(yyvsp[-1].c)+1)*sizeof(char));
       int i;
@@ -8752,13 +8751,13 @@ case 467:
     ;
     break;}
 case 468:
-#line 3869 "Gmsh.y"
+#line 3868 "Gmsh.y"
 {
       yyval.c = yyvsp[-1].c;
     ;
     break;}
 case 469:
-#line 3873 "Gmsh.y"
+#line 3872 "Gmsh.y"
 {
       char tmpstring[1024];
       int i = PrintListOfDouble(yyvsp[-3].c, yyvsp[-1].l, tmpstring);
@@ -8779,7 +8778,7 @@ case 469:
     ;
     break;}
 case 470:
-#line 3892 "Gmsh.y"
+#line 3891 "Gmsh.y"
 { 
       char* (*pStrOpt)(int num, int action, char *value);
       StringXString *pStrCat;
@@ -8797,7 +8796,7 @@ case 470:
     ;
     break;}
 case 471:
-#line 3908 "Gmsh.y"
+#line 3907 "Gmsh.y"
 { 
       char* (*pStrOpt)(int num, int action, char *value);
       StringXString *pStrCat;
@@ -9036,7 +9035,7 @@ yyerrhandle:
     }
   return 1;
 }
-#line 3925 "Gmsh.y"
+#line 3924 "Gmsh.y"
 
 
 void DeleteSymbol(void *a, void *b){
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 8ee8aa4162..f299a3a637 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1,5 +1,5 @@
 %{
-// $Id: Gmsh.y,v 1.167 2004-05-25 04:10:10 geuzaine Exp $
+// $Id: Gmsh.y,v 1.168 2004-05-25 23:16:32 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
@@ -3208,7 +3208,6 @@ Transfinite :
         Surface *s = FindSurface(j, THEM);
 	if(s){
 	  s->Recombine = 1;
-	  s->RecombineAngle = 30.;
         }
       }
       List_Delete($3);
diff --git a/Parser/Gmsh.yy.cpp b/Parser/Gmsh.yy.cpp
index e3b220c889..1afcc7eae9 100644
--- a/Parser/Gmsh.yy.cpp
+++ b/Parser/Gmsh.yy.cpp
@@ -2,7 +2,7 @@
 /* A lexical scanner generated by flex */
 
 /* Scanner skeleton version:
- * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.190 2004-05-25 04:10:10 geuzaine Exp $
+ * $Header: /cvsroot/gmsh/Parser/Gmsh.yy.cpp,v 1.191 2004-05-25 23:16:32 geuzaine Exp $
  */
 
 #define FLEX_SCANNER
@@ -1014,7 +1014,7 @@ char *yytext;
 #line 1 "Gmsh.l"
 #define INITIAL 0
 #line 2 "Gmsh.l"
-// $Id: Gmsh.yy.cpp,v 1.190 2004-05-25 04:10:10 geuzaine Exp $
+// $Id: Gmsh.yy.cpp,v 1.191 2004-05-25 23:16:32 geuzaine Exp $
 //
 // Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 //
diff --git a/Parser/Makefile b/Parser/Makefile
index 127f7b662d..527e0db3b9 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.57 2004-05-12 02:02:30 geuzaine Exp $
+# $Id: Makefile,v 1.58 2004-05-25 23:16:32 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -63,15 +63,15 @@ Gmsh.yy.o: Gmsh.yy.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
   ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h Gmsh.tab.hpp
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h Gmsh.tab.hpp
 Gmsh.tab.o: Gmsh.tab.cpp ../Plugin/PluginManager.h ../Plugin/Plugin.h \
   ../Common/Options.h ../Common/Message.h ../Common/Views.h \
   ../Common/ColorTable.h ../DataStr/List.h ../Parallel/ParUtil.h \
   ../Common/Gmsh.h ../DataStr/Malloc.h ../DataStr/Tree.h ../DataStr/avl.h \
   ../DataStr/Tools.h ../Numeric/Numeric.h ../Common/Context.h \
   ../Geo/Geo.h ../Geo/GeoUtils.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Mesh/Nurbs.h ../Geo/CAD.h ../Graphics/Draw.h ../Mesh/Create.h \
   ../Geo/StepGeomDatabase.h ../Common/Colors.h Parser.h OpenFile.h \
@@ -82,8 +82,8 @@ OpenFile.o: OpenFile.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h \
   ../Common/Context.h Parser.h OpenFile.h ../Common/CommandLine.h \
   ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/STL.h \
-  ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Views.h \
+  ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
+  ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h ../Common/Views.h \
   ../Common/ColorTable.h ../Geo/MinMax.h ../Common/Visibility.h \
   ../Graphics/ReadImg.h ../Common/GmshUI.h ../Graphics/Draw.h \
   ../Fltk/GUI.h ../Fltk/Opengl_Window.h ../Fltk/Colorbar_Window.h
diff --git a/Plugin/Makefile b/Plugin/Makefile
index e999a9df11..deb4f971ef 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.47 2004-05-12 02:02:30 geuzaine Exp $
+# $Id: Makefile,v 1.48 2004-05-25 23:16:32 geuzaine Exp $
 #
 # Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 #
@@ -109,7 +109,7 @@ Triangulate.o: Triangulate.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/avl.h ../DataStr/Tools.h Plugin.h ../Common/Options.h \
   ../Common/Views.h ../Common/ColorTable.h Triangulate.h \
   ../Common/Context.h ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Edge.h \
+  ../Mesh/Element.h ../Mesh/Simplex.h ../Mesh/Face.h ../Mesh/Edge.h \
   ../Geo/ExtrudeParams.h ../Mesh/STL.h ../Mesh/Metric.h ../Mesh/Matrix.h \
   ../Mesh/Utils.h ../Mesh/Create.h
 SphericalRaise.o: SphericalRaise.cpp Plugin.h ../Common/Options.h \
diff --git a/doc/VERSIONS b/doc/VERSIONS
index d87247ad02..abd918e620 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,16 +1,17 @@
-$Id: VERSIONS,v 1.210 2004-05-25 04:10:10 geuzaine Exp $
-
-New in 1.53: various background mesh fixes and enhancements; new
-Plugin(Evaluate) to evaluate arbitrary expressions on post-processing
-views; generalized Plugin(Extract) to handle any combination of
-components; generalized "Coherence" to handle transfinite
-surface/volume attributes; plugin options can now be set in the option
-file (like all other options); added "undo" capability during geometry
-creation; rewrote the contour guessing routines so that entities can
-be selected in an arbitrary order; Mac users can now double click on
-geo/msh/pos files in the Finder to launch Gmsh; removed support for
-fltk 1.0; rewrote most of the code related to quadrangles; many code
-cleanups;
+$Id: VERSIONS,v 1.211 2004-05-25 23:16:32 geuzaine Exp $
+
+New in 1.53: completed support for second order elements (lines,
+triangles, quadrangles, tetrahedra, hexahedra, prisms and pyramids);
+various background mesh fixes and enhancements; new Plugin(Evaluate)
+to evaluate arbitrary expressions on post-processing views;
+generalized Plugin(Extract) to handle any combination of components;
+generalized "Coherence" to handle transfinite surface/volume
+attributes; plugin options can now be set in the option file (like all
+other options); added "undo" capability during geometry creation;
+rewrote the contour guessing routines so that entities can be selected
+in an arbitrary order; Mac users can now double click on geo/msh/pos
+files in the Finder to launch Gmsh; removed support for fltk 1.0;
+rewrote most of the code related to quadrangles; many code cleanups;
 
 New in 1.52: new raster ("bitmap") PostScript/EPS/PDF output formats;
 new Plugin(Extract) to extract a given component from a
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index ee474f9f5f..702f625885 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -1,5 +1,5 @@
 \input texinfo.tex @c -*-texinfo-*-
-@c $Id: gmsh.texi,v 1.111 2004-05-25 04:20:34 geuzaine Exp $
+@c $Id: gmsh.texi,v 1.112 2004-05-25 23:16:33 geuzaine Exp $
 @c
 @c Copyright (C) 1997-2004 C. Geuzaine, J.-F. Remacle
 @c
@@ -2737,9 +2737,9 @@ is the list of the @var{number-of-nodes} node numbers of the @var{n}-th
 element (separated by white space, without commas). The ordering of the
 nodes is given in @ref{Gmsh node ordering}; for second order elements, the
 first order nodes are given first, followed by the nodes associated with the
-edges, followed by the nodes associated with the faces (if any). The
-ordering of these additional nodes follows the ordering of the edges/faces
-given in @ref{Gmsh node ordering}.
+edges, followed by the nodes associated with the quadrangular faces (if
+any). The ordering of these additional nodes follows the ordering of the
+edges/faces given in @ref{Gmsh node ordering}.
 @end table
 
 @c .........................................................................
@@ -2849,9 +2849,9 @@ is the list of the node numbers of the @var{n}-th element (separated by
 white space, without commas). The ordering of the nodes is given in
 @ref{Gmsh node ordering}; for second order elements, the first order nodes
 are given first, followed by the nodes associated with the edges, followed
-by the nodes associated with the faces (if any). The ordering of these
-additional nodes follows the ordering of the edges/faces given in @ref{Gmsh
-node ordering}.
+by the nodes associated with the quadrangular faces (if any). The ordering
+of these additional nodes follows the ordering of the edges/faces given in
+@ref{Gmsh node ordering}.
 @end table
 
 @c -------------------------------------------------------------------------
@@ -3172,102 +3172,102 @@ defined as follows:
 @example
 @group
 Point:
-        v                            
-        |                            
-        |                            
-   -----1-----u  
-        |                            
-        |                            
+        v
+        |
+        |
+   -----1-----u
+        |
+        |
 @end group
 @end example
 
 @example
 @group
 Line:
-                  edge 1: nodes 1 -> 2  
-        v                            
-        |                            
-        |                            
-   --1-----2--u  
-        |                            
-        |                            
+                  edge 1: nodes 1 -> 2
+        v
+        |
+        |
+   --1-----2--u
+        |
+        |
 @end group
 @end example
 
 @example
 @group
 Triangle:
-                  edge 1: nodes 1 -> 2         
-   v                   2:       2 -> 3                  
-   |                   3:       3 -> 1                  
-   |                           
-   3              
-   |\        
-   | \       
-   |__\___u                    
-   1   2                       
+                  edge 1: nodes 1 -> 2
+   v                   2:       2 -> 3
+   |                   3:       3 -> 1
+   |
+   3
+   |\
+   | \
+   |__\___u
+   1   2
 @end group
 @end example
 
 @example
 @group
 Quadrangle:
-                  edge 1: nodes 1 -> 2 
-        v              2:       2 -> 3                    
-        |              3:       3 -> 4                    
-     4--|--3           4:       4 -> 1 
-     |  |  |     
-   -----------u   
-     |  |  |     
-     1--|--2      
-        |         
+                  edge 1: nodes 1 -> 2   quadface 1: nodes 1 2 3 4
+        v              2:       2 -> 3
+        |              3:       3 -> 4
+     4--|--3           4:       4 -> 1
+     |  |  |
+   -----------u
+     |  |  |
+     1--|--2
+        |
 @end group
 @end example
 
 @example
 @group
 Tetrahedron:
-                  edge 1: nodes 1 -> 2                
-   v                   2:       2 -> 3                                            
-   |                   3:       3 -> 1                
-   |                   4:       4 -> 1                
-   |                   5:       4 -> 3                
-   3                   6:       4 -> 2                
-   |\                                                 
-   | \            
-   |__\2_____u    
-   1\ /           
-     \4           
-      \          
-       w         
+                  edge 1: nodes 1 -> 2
+   v                   2:       2 -> 3
+   |                   3:       3 -> 1
+   |                   4:       4 -> 1
+   |                   5:       4 -> 3
+   3                   6:       4 -> 2
+   |\                
+   | \
+   |__\2_____u
+   1\ /
+     \4
+      \
+       w
 @end group
 @end example
 
 @example
 @group
 Hexahedron:
-                  edge 1: nodes 1 -> 2
-        v              2:       1 -> 4
-        |              3:       1 -> 5
-        |              4:       2 -> 3
-   4----|--3           5:       2 -> 6
-   |\   |  |\          6:       3 -> 4
+                  edge 1: nodes 1 -> 2   quadface 1: nodes 1 2 3 4
+        v              2:       1 -> 4            2:       1 2 5 6
+        |              3:       1 -> 5            3:       1 4 5 8
+        |              4:       2 -> 3            4:       2 3 6 7
+   4----|--3           5:       2 -> 6            5:       3 4 7 8
+   |\   |  |\          6:       3 -> 4            6:       5 6 7 8
    | 8-------7         7:       3 -> 7
    | |   ----|---u     8:       4 -> 8
    1-|---\-2 |         9:       5 -> 6
     \|    \ \|        10:       5 -> 8
      5-----\-6        11:       6 -> 7
             \         12:       7 -> 8
-             w         
+             w
 @end group
 @end example
 
 @example
 @group
 Prism:
-                  edge 1: nodes 1 -> 2
-      v                2:       1 -> 3
-    3 |                3:       1 -> 4
+                  edge 1: nodes 1 -> 2   quadface 1: nodes 1 2 4 5
+      v                2:       1 -> 3            2:       1 3 4 6
+    3 |                3:       1 -> 4            3:       2 3 5 6
     |\|                4:       2 -> 3
     | |                5:       2 -> 5
     1_|2               6:       3 -> 6
@@ -3275,16 +3275,16 @@ Prism:
       |_|_\___u        8:       4 -> 6
        \|  \           9:       5 -> 6
         4 __5
-         \       
-          \      
-           w     
+         \
+          \
+           w
 @end group
 @end example
 
 @example
 @group
 Pyramid:
-                  edge 1: nodes 1 -> 2
+                  edge 1: nodes 1 -> 2   quadface 1: nodes 1 2 3 4
         v              2        1 -> 4
         |              3        1 -> 5
         |              4        2 -> 3
@@ -3292,10 +3292,10 @@ Pyramid:
     | \ |  /|          6        3 -> 4
     |  \ -/-|---u      7        3 -> 5
     |  / 5\ |          8        4 -> 5
-    1/----\-2              
-           \      
-            \     
-             w    
+    1/----\-2
+           \
+            \
+             w
 @end group
 @end example
 
-- 
GitLab