diff --git a/Common/Context.h b/Common/Context.h
index cc0edf3569305e85cfa2c7ab37f6a1b08d86603e..d4fc3f081f3f4dcb1213cf958cff493fb0ec2059 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -170,7 +170,8 @@ public :
     double scaling_factor, lc_factor, rand_factor;
     int dual, interactive;
     int light, light_two_side, light_lines;
-    int format, nbPartitions, nb_smoothing, algo2d, algo3d, order, algo_recombine;
+    int format, nbPartitions, nb_smoothing, algo2d, algo3d, algo_recombine;
+    int order, second_order_linear, second_order_incomplete;
     int point_insertion, speed_max, min_circ_points;
     int bgmesh_type, constrained_bgmesh;
     int initial_only;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index b7fda784a228ecd5bc2215997805ee750c890d70..62d046b4590d66a7e0b02cdb22ed3019a892f256 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -943,6 +943,8 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "Optimize" , opt_mesh_optimize , 0. , 
     "Optimize the mesh using Netgen to improve the quality of tetrahedral elements" },
 
+  { F|O, "Partitioning" , opt_mesh_nb_partitions , 1. ,
+    "Number of partitions applied to the final mesh" },
   { F|O, "Points" , opt_mesh_points , 0. , 
     "Display mesh vertices (nodes)?" },
   { F|O, "PointInsertion" , opt_mesh_point_insertion, CENTER_CIRCCIRC ,
@@ -984,10 +986,12 @@ StringXNumber MeshOptions_Number[] = {
     "Ignore Physical definitions and save all elements" },
   { F|O, "ScalingFactor" , opt_mesh_scaling_factor , 1.0 ,
     "Global scaling factor applied to the saved mesh" },
+  { F|O, "SecondOrderIncomplete" , opt_mesh_second_order_incomplete , 1. ,
+    "Create incomplete second order elements? (8-node quads, 20-node hexas, etc.)" },
+  { F|O, "SecondOrderLinear" , opt_mesh_second_order_linear , 1. ,
+    "Should second order vertices simply be created by linear interpolation?" },
   { F|O, "Smoothing" , opt_mesh_nb_smoothing , 1. ,
     "Number of smoothing steps applied to the final mesh" },
-  { F|O, "Partitioning" , opt_mesh_nb_partitions , 1. ,
-    "Number of partitions applied to the final mesh" },
   { F|O, "SmoothNormals" , opt_mesh_smooth_normals , 0. , 
     "Smooth the mesh normals?" },
   { F|O, "SpeedMax" , opt_mesh_speed_max , 0. ,
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 41f7f4bc05dbdcc90f28ec14e43b78ccfcf9b5ea..43419128ef053c3b2a2c12704bf1dd25fe0d2171 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -31,21 +31,25 @@
 #define FORMAT_BDF           31
 
 // Element types in .msh file format
-#define LGN1 1
-#define TRI1 2
-#define QUA1 3
-#define TET1 4
-#define HEX1 5
-#define PRI1 6
-#define PYR1 7
-#define LGN2 8
-#define TRI2 9
-#define QUA2 10
-#define TET2 11
-#define HEX2 12
-#define PRI2 13
-#define PYR2 14
-#define PNT  15
+#define LIN_2  1
+#define TRI_3  2
+#define QUA_4  3
+#define TET_4  4
+#define HEX_8  5
+#define PRI_6  6
+#define PYR_5  7
+#define LIN_3  8
+#define TRI_6  9
+#define QUA_9  10
+#define TET_10 11
+#define HEX_27 12
+#define PRI_18 13
+#define PYR_14 14
+#define PNT_1  15
+#define QUA_8  16
+#define HEX_20 17
+#define PRI_15 18
+#define PYR_13 19
 
 // Geometric entities
 #define ENT_NONE     0
diff --git a/Common/Options.cpp b/Common/Options.cpp
index d6272aa239f9029a8083c44be123c6e7c47cbdcf..f38860f7114a1295e7b42c81363c5abaacca5965 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.309 2006-09-03 07:44:09 geuzaine Exp $
+// $Id: Options.cpp,v 1.310 2006-09-08 02:39:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -4829,6 +4829,20 @@ double opt_mesh_order(OPT_ARGS_NUM)
   return CTX.mesh.order;
 }
 
+double opt_mesh_second_order_linear(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.mesh.second_order_linear = (int)val;
+  return CTX.mesh.second_order_linear;
+}
+
+double opt_mesh_second_order_incomplete(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.mesh.second_order_incomplete = (int)val;
+  return CTX.mesh.second_order_incomplete;
+}
+
 double opt_mesh_dual(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET) {
diff --git a/Common/Options.h b/Common/Options.h
index 5ae072b19e24b293a75697e934fe391676bcaf61..d070e04baa05dfa932a732316e1f88223780b9f8 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -466,6 +466,8 @@ double opt_mesh_speed_max(OPT_ARGS_NUM);
 double opt_mesh_min_circ_points(OPT_ARGS_NUM);
 double opt_mesh_constrained_bgmesh(OPT_ARGS_NUM);
 double opt_mesh_order(OPT_ARGS_NUM);
+double opt_mesh_second_order_linear(OPT_ARGS_NUM);
+double opt_mesh_second_order_incomplete(OPT_ARGS_NUM);
 double opt_mesh_dual(OPT_ARGS_NUM);
 double opt_mesh_interactive(OPT_ARGS_NUM);
 double opt_mesh_initial_only(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 55463e2c1bbd9dfa8e187df909854aece3455bdb..bff3115816ec0a3eb42dbb27ae2131ba520f7174 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.460 2006-09-03 07:44:10 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.461 2006-09-08 02:39:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -3505,7 +3505,7 @@ void mesh_reparam_cb(CALLBACK_ARGS)
 void mesh_degree_cb(CALLBACK_ARGS)
 {
   if((long)data == 2)
-    Degre2();
+    Degre2(CTX.mesh.second_order_linear, CTX.mesh.second_order_incomplete);
   else
     Degre1();
   CTX.mesh.changed = ENT_LINE | ENT_SURFACE | ENT_VOLUME;
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index a540a04f25263b59044d9b03fe69980c394a7b5b..f38ba9323fd887b1c627f52ce7dce2546d5902cf 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.547 2006-08-31 23:01:11 geuzaine Exp $
+// $Id: GUI.cpp,v 1.548 2006-09-08 02:39:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -3640,7 +3640,7 @@ void GUI::create_statistics_window()
   }
 
   int width = 26 * fontsize;
-  int height = 5 * WB + 18 * BH;
+  int height = 5 * WB + 17 * BH;
 
   stat_window = new Dialog_Window(width, height, "Statistics");
   stat_window->box(GMSH_WINDOW_BOX);
@@ -3660,27 +3660,26 @@ void GUI::create_statistics_window()
       stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 1 * BH, IW, BH, "Nodes on Lines");
       stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 2 * BH, IW, BH, "Nodes on surfaces");
       stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 3 * BH, IW, BH, "Nodes in volumes");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 4 * BH, IW, BH, "Second order nodes");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 5 * BH, IW, BH, "Triangles");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 6 * BH, IW, BH, "Quadrangles");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 7 * BH, IW, BH, "Tetrahedra");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 8 * BH, IW, BH, "Hexahedra");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 9 * BH, IW, BH, "Prisms");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Pyramids");
-
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Time for 1D mesh");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 12 * BH, IW, BH, "Time for 2D mesh");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "Time for 3D mesh");
-
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "Gamma");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Eta");
-      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 16 * BH, IW, BH, "Rho");
-
-      stat_butt[0] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 14 * BH, BB, BH, "Graph");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 4 * BH, IW, BH, "Triangles");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 5 * BH, IW, BH, "Quadrangles");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 6 * BH, IW, BH, "Tetrahedra");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 7 * BH, IW, BH, "Hexahedra");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 8 * BH, IW, BH, "Prisms");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 9 * BH, IW, BH, "Pyramids");
+
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 10 * BH, IW, BH, "Time for 1D mesh");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 11 * BH, IW, BH, "Time for 2D mesh");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 12 * BH, IW, BH, "Time for 3D mesh");
+
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 13 * BH, IW, BH, "Gamma");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 14 * BH, IW, BH, "Eta");
+      stat_value[num++] = new Fl_Output(2 * WB, 2 * WB + 15 * BH, IW, BH, "Rho");
+
+      stat_butt[0] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 13 * BH, BB, BH, "Graph");
       stat_butt[0]->callback(statistics_histogram_cb, (void *)"Gamma");
-      stat_butt[1] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 15 * BH, BB, BH, "Graph");
+      stat_butt[1] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 14 * BH, BB, BH, "Graph");
       stat_butt[1]->callback(statistics_histogram_cb, (void *)"Eta");
-      stat_butt[2] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 16 * BH, BB, BH, "Graph");
+      stat_butt[2] = new Fl_Button(width - BB - 5 * WB, 2 * WB + 15 * BH, BB, BH, "Graph");
       stat_butt[2]->callback(statistics_histogram_cb, (void *)"Rho");
 
       g[1]->end();
@@ -3716,7 +3715,6 @@ void GUI::create_statistics_window()
     o->callback(cancel_cb, (void *)stat_window);
   }
 
-
   stat_window->position(CTX.stat_position[0], CTX.stat_position[1]);
   stat_window->end();
 }
@@ -3745,7 +3743,6 @@ void GUI::set_statistics(bool compute_quality)
   sprintf(label[num], "%g", s[4]); stat_value[num]->value(label[num]); num++;
   sprintf(label[num], "%g", s[5]); stat_value[num]->value(label[num]); num++;
   sprintf(label[num], "%g", s[6]); stat_value[num]->value(label[num]); num++;
-  sprintf(label[num], "%g", s[16]); stat_value[num]->value(label[num]); num++;
   sprintf(label[num], "%g", s[7]); stat_value[num]->value(label[num]); num++;
   sprintf(label[num], "%g", s[8]); stat_value[num]->value(label[num]); num++;
   sprintf(label[num], "%g", s[9]); stat_value[num]->value(label[num]); num++;
diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index a49769ed0d61407a45498aa0b73b134680a10a94..08299a4bd0cf366597122045dadf6c6e30db76e3 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.44 2006-09-07 19:45:15 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.45 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -194,21 +194,25 @@ static bool checkVertexIndex(int index, std::vector<MVertex*> &vec)
 static int getNumVerticesForElementTypeMSH(int type)
 {
   switch (type) {
-  case PNT  : return 1;
-  case LGN1 : return 2;
-  case LGN2 : return 2 + 1;
-  case TRI1 : return 3;
-  case TRI2 : return 3 + 3;
-  case QUA1 : return 4;
-  case QUA2 : return 4 + 4 + 1;
-  case TET1 : return 4;
-  case TET2 : return 4 + 6;
-  case HEX1 : return 8;
-  case HEX2 : return 8 + 12 + 6 + 1;
-  case PRI1 : return 6;
-  case PRI2 : return 6 + 9 + 3;
-  case PYR1 : return 5;
-  case PYR2 : return 5 + 8 + 1;
+  case PNT_1  : return 1;
+  case LIN_2  : return 2;
+  case LIN_3  : return 2 + 1;
+  case TRI_3  : return 3;
+  case TRI_6  : return 3 + 3;
+  case QUA_4  : return 4;
+  case QUA_8  : return 4 + 4;
+  case QUA_9  : return 4 + 4 + 1;
+  case TET_4  : return 4;
+  case TET_10 : return 4 + 6;
+  case HEX_8  : return 8;
+  case HEX_20 : return 8 + 12;
+  case HEX_27 : return 8 + 12 + 6 + 1;
+  case PRI_6  : return 6;
+  case PRI_15 : return 6 + 9;
+  case PRI_18 : return 6 + 9 + 3;
+  case PYR_5  : return 5;
+  case PYR_13 : return 5 + 8;
+  case PYR_14 : return 5 + 8 + 1;
   default: 
     Msg(GERROR, "Unknown type of element %d", type);
     return 0;
@@ -225,60 +229,74 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
   int dim = 0;
 
   switch (type) {
-  case PNT:
+  case PNT_1:
     points[elementary].push_back(v[n[0]]);
     dim = 0;
     break;
-  case LGN1:
+  case LIN_2:
     elements[0][elementary].push_back
       (new MLine(v[n[0]], v[n[1]], num, partition));
     dim = 1;
     break;
-  case LGN2:
+  case LIN_3:
     elements[0][elementary].push_back
       (new MLine3(v[n[0]], v[n[1]], v[n[2]], num, partition));
     dim = 1;
     break;
-  case TRI1:
+  case TRI_3:
     elements[1][elementary].push_back
       (new MTriangle(v[n[0]], v[n[1]], v[n[2]], num, partition));
     dim = 2;
     break;
-  case TRI2:
+  case TRI_6:
     elements[1][elementary].push_back
       (new MTriangle6(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 		      num, partition));
     dim = 2;
     break;
-  case QUA1:
+  case QUA_4:
     elements[2][elementary].push_back
       (new MQuadrangle(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
     dim = 2;
     break;
-  case QUA2:
+  case QUA_8:
+    elements[2][elementary].push_back
+      (new MQuadrangle8(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+			v[n[6]], v[n[7]], num, partition));
+    dim = 2;
+    break;
+  case QUA_9:
     elements[2][elementary].push_back
       (new MQuadrangle9(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 			v[n[6]], v[n[7]], v[n[8]], num, partition));
     dim = 2;
     break;
-  case TET1:
+  case TET_4:
     elements[3][elementary].push_back
       (new MTetrahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
     dim = 3; 
     break;
-  case TET2:
+  case TET_10:
     elements[3][elementary].push_back
       (new MTetrahedron10(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 			  v[n[6]], v[n[7]], v[n[8]], v[n[9]], num, partition));
     dim = 3; 
     break;
-  case HEX1:
+  case HEX_8:
     elements[4][elementary].push_back
       (new MHexahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 		       v[n[6]], v[n[7]], num, partition));
     dim = 3; 
     break;
-  case HEX2:
+  case HEX_20:
+    elements[4][elementary].push_back
+      (new MHexahedron20(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+			 v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+			 v[n[12]], v[n[13]], v[n[14]], v[n[15]], v[n[16]], v[n[17]], 
+			 v[n[18]], v[n[19]], num, partition));
+    dim = 3; 
+    break;
+  case HEX_27:
     elements[4][elementary].push_back
       (new MHexahedron27(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 			 v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
@@ -287,13 +305,20 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
 			 v[n[24]], v[n[25]], v[n[26]], num, partition));
     dim = 3; 
     break;
-  case PRI1: 
+  case PRI_6: 
     elements[5][elementary].push_back
       (new MPrism(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 		  num, partition));
     dim = 3; 
     break;
-  case PRI2: 
+  case PRI_15: 
+    elements[5][elementary].push_back
+      (new MPrism15(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		    v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+		    v[n[12]], v[n[13]], v[n[14]], num, partition));
+    dim = 3; 
+    break;
+  case PRI_18: 
     elements[5][elementary].push_back
       (new MPrism18(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 		    v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
@@ -301,12 +326,19 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
 		    num, partition));
     dim = 3; 
     break;
-  case PYR1: 
+  case PYR_5: 
     elements[6][elementary].push_back
       (new MPyramid(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], num, partition));
     dim = 3; 
     break;
-  case PYR2: 
+  case PYR_13: 
+    elements[6][elementary].push_back
+      (new MPyramid13(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		      v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+		      v[n[12]], num, partition));
+    dim = 3; 
+    break;
+  case PYR_14: 
     elements[6][elementary].push_back
       (new MPyramid14(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
 		      v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
@@ -558,10 +590,25 @@ int GModel::readMSH(const std::string &name)
   return 1;
 }
 
-static void writeTagMSH(FILE *fp, int type, int num, int numTags)
+static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elements,
+				  int t1, int t2=0, int t3=0)
 {
-  int data[3] = {type, num, numTags};
-  fwrite(data, sizeof(int), 3, fp);
+  if(!binary) return;
+
+  int numTags = 3;
+  int data[3];
+  if(elements.count(t1)){
+    data[0] = t1;  data[1] = elements[t1];  data[2] = numTags;
+    fwrite(data, sizeof(int), 3, fp);
+  }
+  else if(t2 && elements.count(t2)){
+    data[0] = t2;  data[1] = elements[t2];  data[2] = numTags;
+    fwrite(data, sizeof(int), 3, fp);
+  }
+  else if(t3 && elements.count(t3)){
+    data[0] = t3;  data[1] = elements[t3];  data[2] = numTags;
+    fwrite(data, sizeof(int), 3, fp);
+  }
 }
 
 template<class T>
@@ -603,7 +650,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   for(viter it = firstVertex(); it != lastVertex(); ++it){
     int p = (saveAll ? 1 : (*it)->physicals.size());
     int n = p * (*it)->mesh_vertices.size();
-    if(n) elements[PNT] += n;
+    if(n) elements[PNT_1] += n;
   }
   for(eiter it = firstEdge(); it != lastEdge(); ++it){
     int p = (saveAll ? 1 : (*it)->physicals.size());
@@ -674,51 +721,37 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   }
 
   fprintf(fp, "%d\n", numElements);
-  int num = 0, numTags = 3;
+  int num = 0;
 
-  if(binary && elements.count(PNT)) writeTagMSH(fp, PNT, elements[PNT], numTags);
+  writeElementHeaderMSH(binary, fp, elements, PNT_1);
   for(viter it = firstVertex(); it != lastVertex(); ++it)
     writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(LGN1)) writeTagMSH(fp, LGN1, elements[LGN1], numTags);
-  else if(binary && elements.count(LGN2)) writeTagMSH(fp, LGN2, elements[LGN2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, LIN_2, LIN_3);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
     writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(TRI1)) writeTagMSH(fp, TRI1, elements[TRI1], numTags);
-  else if(binary && elements.count(TRI2)) writeTagMSH(fp, TRI2, elements[TRI2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, TRI_3, TRI_6);
   for(fiter it = firstFace(); it != lastFace(); ++it)
     writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(QUA1)) writeTagMSH(fp, QUA1, elements[QUA1], numTags);
-  else if(binary && elements.count(QUA2)) writeTagMSH(fp, QUA2, elements[QUA2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, QUA_4, QUA_9, QUA_8);
   for(fiter it = firstFace(); it != lastFace(); ++it)
     writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(TET1)) writeTagMSH(fp, TET1, elements[TET1], numTags);
-  else if(binary && elements.count(TET2)) writeTagMSH(fp, TET2, elements[TET2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, TET_4, TET_10);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(HEX1)) writeTagMSH(fp, HEX1, elements[HEX1], numTags);
-  else if(binary && elements.count(HEX2)) writeTagMSH(fp, HEX2, elements[HEX2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, HEX_8, HEX_27, HEX_20);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(PRI1)) writeTagMSH(fp, PRI1, elements[PRI1], numTags);
-  else if(binary && elements.count(PRI2)) writeTagMSH(fp, PRI2, elements[PRI2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, PRI_6, PRI_18, PRI_15);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-
-  if(binary && elements.count(PYR1)) writeTagMSH(fp, PYR1, elements[PYR1], numTags);
-  else if(binary && elements.count(PYR2)) writeTagMSH(fp, PYR2, elements[PYR2], numTags);
+  writeElementHeaderMSH(binary, fp, elements, PYR_5, PYR_14, PYR_13);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
@@ -1217,7 +1250,11 @@ int GModel::readUNV(const std::string &name)
 	    dim = 2;
 	    break;
 	  case 45: case 55: case 65: case 75: case 85: case 95:
-	    Msg(WARNING, "Quad8 is not implemented yet");
+	    elements[2][elementary].push_back
+	      (new MQuadrangle8(vertices[n[0]], vertices[n[2]], vertices[n[4]], 
+				vertices[n[6]], vertices[n[1]], vertices[n[3]], 
+				vertices[n[5]], vertices[n[7]], num));
+	    dim = 2;
 	    break;
 	  case 111:
 	    elements[3][elementary].push_back
@@ -1241,7 +1278,15 @@ int GModel::readUNV(const std::string &name)
 	    dim = 3; 
 	    break;
 	  case 105: case 116:
-	    Msg(WARNING, "Hexa20 is not implemented yet");
+	    elements[4][elementary].push_back
+	      (new MHexahedron20(vertices[n[0]], vertices[n[2]], vertices[n[4]],  
+				 vertices[n[6]], vertices[n[12]], vertices[n[14]], 
+				 vertices[n[16]], vertices[n[18]], vertices[n[1]],  
+				 vertices[n[7]], vertices[n[8]], vertices[n[3]],  
+				 vertices[n[9]], vertices[n[5]], vertices[n[10]], 
+				 vertices[n[11]], vertices[n[13]], vertices[n[19]], 
+				 vertices[n[15]], vertices[n[17]], num));
+	    dim = 3; 
 	    break;
 	  case 101: case 112:
 	    elements[5][elementary].push_back
@@ -1250,7 +1295,13 @@ int GModel::readUNV(const std::string &name)
 	    dim = 3; 
 	    break;
 	  case 102: case 113:
-	    Msg(WARNING, "Prism15 is not implemented yet");
+	    elements[5][elementary].push_back
+	      (new MPrism15(vertices[n[0]], vertices[n[2]], vertices[n[4]], 
+			    vertices[n[9]], vertices[n[11]], vertices[n[13]], 
+			    vertices[n[1]], vertices[n[5]], vertices[n[6]],  
+			    vertices[n[3]], vertices[n[7]], vertices[n[8]],  
+			    vertices[n[10]], vertices[n[14]], vertices[n[12]], num));
+	    dim = 3;
 	    break;
 	  }
   
@@ -1677,6 +1728,7 @@ int GModel::readBDF(const std::string &name)
 	  elements[1][region].push_back(new MTriangle(n, num));
       }
       else if(!strncmp(buffer, "CTRIA6", 6)){
+	// not sure about the node ordering
 	if(readElementBDF(fp, buffer, 6, 6, &num, &region, n, vertices))
 	  elements[1][region].push_back(new MTriangle6(n, num));
       }
@@ -1684,6 +1736,11 @@ int GModel::readBDF(const std::string &name)
 	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertices))
 	  elements[2][region].push_back(new MQuadrangle(n, num));
       }
+      else if(!strncmp(buffer, "CQUAD8", 6)){
+	// not sure about the node ordering
+	if(readElementBDF(fp, buffer, 6, 8, &num, &region, n, vertices))
+	  elements[2][region].push_back(new MQuadrangle8(n, num));
+      }
       else if(!strncmp(buffer, "CTETRA", 6)){
 	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertices))
 	  elements[3][region].push_back(new MTetrahedron(n, num));
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 08ce1e8bff86e6eaf415a5d811b5ed9410d4f684..9ab1d5b82075f851333cc88fb05cc6905041ef2f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.16 2006-09-07 19:45:15 geuzaine Exp $
+// $Id: MElement.cpp,v 1.17 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -163,33 +163,35 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
 
 void MElement::writePOS(FILE *fp, double scalingFactor, int elementary)
 {
-  int n = getNumVertices();
-  double gamma = gammaShapeMeasure();
-  double eta = etaShapeMeasure();
-  double rho = rhoShapeMeasure();
-
-  fprintf(fp, "%s(", getStringForPOS());
-  for(int i = 0; i < n; i++){
-    if(i) fprintf(fp, ",");
-    fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor, 
-	    getVertex(i)->y() * scalingFactor, getVertex(i)->z() * scalingFactor);
-  }
-  fprintf(fp, "){");
-  for(int i = 0; i < n; i++)
-    fprintf(fp, "%d,", elementary);
-  for(int i = 0; i < n; i++)
-    fprintf(fp, "%d,", getNum());
-  for(int i = 0; i < n; i++)
-    fprintf(fp, "%g,", gamma);
-  for(int i = 0; i < n; i++)
-    fprintf(fp, "%g,", eta);
-  for(int i = 0; i < n; i++){
-    if(i == n - 1)
-      fprintf(fp, "%g", rho);
-    else
-      fprintf(fp, "%g,", rho);
+  char *str = getStringForPOS();
+  if(str){
+    int n = getNumVertices();
+    double gamma = gammaShapeMeasure();
+    double eta = etaShapeMeasure();
+    double rho = rhoShapeMeasure();
+    fprintf(fp, "%s(", str);
+    for(int i = 0; i < n; i++){
+      if(i) fprintf(fp, ",");
+      fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor, 
+	      getVertex(i)->y() * scalingFactor, getVertex(i)->z() * scalingFactor);
+    }
+    fprintf(fp, "){");
+    for(int i = 0; i < n; i++)
+      fprintf(fp, "%d,", elementary);
+    for(int i = 0; i < n; i++)
+      fprintf(fp, "%d,", getNum());
+    for(int i = 0; i < n; i++)
+      fprintf(fp, "%g,", gamma);
+    for(int i = 0; i < n; i++)
+      fprintf(fp, "%g,", eta);
+    for(int i = 0; i < n; i++){
+      if(i == n - 1)
+	fprintf(fp, "%g", rho);
+      else
+	fprintf(fp, "%g,", rho);
+    }
+    fprintf(fp, "};\n");
   }
-  fprintf(fp, "};\n");
 }
 
 void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index e24ef8e4a8a9ffc16a579d92b5bafa7428f28ce9..d094b4e9fea73636171c69f2cc4b229a98a571b7 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -141,10 +141,12 @@ class MElement
   virtual void writeUNV(FILE *fp, int num=0, int elementary=1, int physical=1);
   virtual void writeMESH(FILE *fp, int elementary=1);
   virtual void writeBDF(FILE *fp, int format=0, int elementary=1);
-  virtual int getTypeForMSH() = 0;
-  virtual int getTypeForUNV() = 0;
-  virtual char *getStringForPOS() = 0;
-  virtual char *getStringForBDF() = 0;
+
+  // default value means that the particular element type is not available
+  virtual int getTypeForMSH(){ return 0; }
+  virtual int getTypeForUNV(){ return 0; }
+  virtual char *getStringForPOS(){ return 0; }
+  virtual char *getStringForBDF(){ return 0; }
 };
 
 class MLine : public MElement {
@@ -169,7 +171,7 @@ class MLine : public MElement {
   virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
   virtual int getNumFaces(){ return 0; }
   virtual MFace getFace(int num){ throw; }
-  virtual int getTypeForMSH(){ return LGN1; }
+  virtual int getTypeForMSH(){ return LIN_2; }
   virtual int getTypeForUNV(){ return 21; } // linear beam
   virtual char *getStringForPOS(){ return "SL"; }
   virtual char *getStringForBDF(){ return "CBAR"; }
@@ -209,10 +211,9 @@ class MLine3 : public MLine {
     int i1 = edges_lin2[num][1];
     return MEdge(i0 < 2 ? _v[i0] : _vs[i0 - 2], i1 < 2 ? _v[i1] : _vs[i1 - 2]);
   }
-  virtual int getTypeForMSH(){ return LGN2; }
+  virtual int getTypeForMSH(){ return LIN_3; }
   virtual int getTypeForUNV(){ return 24; } // parabolic beam
   virtual char *getStringForPOS(){ return "SL2"; }
-  virtual char *getStringForBDF(){ return 0; } // not available
 };
 
 class MTriangle : public MElement {
@@ -249,7 +250,7 @@ class MTriangle : public MElement {
   { 
     return MFace(_v[0], _v[1], _v[2]); 
   }
-  virtual int getTypeForMSH(){ return TRI1; }
+  virtual int getTypeForMSH(){ return TRI_3; }
   virtual int getTypeForUNV(){ return 91; } // thin shell linear triangle
   virtual char *getStringForPOS(){ return "ST"; }
   virtual char *getStringForBDF(){ return "CTRIA3"; }
@@ -308,7 +309,7 @@ class MTriangle6 : public MTriangle {
 		 i1 < 3 ? _v[i1] : _vs[i1 - 3],
 		 i2 < 3 ? _v[i2] : _vs[i2 - 3]);
   }
-  virtual int getTypeForMSH(){ return TRI2; }
+  virtual int getTypeForMSH(){ return TRI_6; }
   virtual int getTypeForUNV(){ return 92; } // thin shell parabolic triangle
   virtual char *getStringForPOS(){ return "ST2"; }
   virtual char *getStringForBDF(){ return "CTRIA6"; }
@@ -339,12 +340,43 @@ class MQuadrangle : public MElement {
   }
   virtual int getNumFaces(){ return 1; }
   virtual MFace getFace(int num){ return MFace(_v[0], _v[1], _v[2], _v[3]); }
-  virtual int getTypeForMSH(){ return QUA1; }
+  virtual int getTypeForMSH(){ return QUA_4; }
   virtual int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
   virtual char *getStringForPOS(){ return "SQ"; }
   virtual char *getStringForBDF(){ return "CQUAD4"; }
 };
 
+class MQuadrangle8 : public MQuadrangle {
+ protected:
+  MVertex *_vs[4];
+ public :
+  MQuadrangle8(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, 
+	       MVertex *v5, MVertex *v6, MVertex *v7, int num=0, int part=0) 
+    : MQuadrangle(v0, v1, v2, v3, num, part)
+  {
+    _vs[0] = v4; _vs[1] = v5; _vs[2] = v6; _vs[3] = v7;
+  }
+  MQuadrangle8(std::vector<MVertex*> &v, int num=0, int part=0) 
+    : MQuadrangle(v, num, part)
+  {
+    for(int i = 0; i < 4; i++) _vs[i] = v[4 + i];
+  }
+  ~MQuadrangle8(){}
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 8; }
+  virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual MVertex *getVertexUNV(int num)
+  {
+    static const int map[8] = {0, 4, 1, 5, 2, 6, 3, 7};
+    return getVertex(map[num]); 
+  }
+  virtual int getNumEdgeVertices(){ return 4; }
+  // TODO: edgeRep, faceRep
+  virtual int getTypeForMSH(){ return QUA_8; }
+  virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
+  virtual char *getStringForBDF(){ return "CQUAD8"; }
+};
+
 class MQuadrangle9 : public MQuadrangle {
  protected:
   MVertex *_vs[5];
@@ -367,16 +399,10 @@ class MQuadrangle9 : public MQuadrangle {
   virtual int getNumEdgeVertices(){ return 4; }
   virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  virtual int getTypeForMSH(){ return QUA2; }
-  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual int getTypeForMSH(){ return QUA_9; }
   virtual char *getStringForPOS(){ return "SQ2"; }
-  virtual char *getStringForBDF(){ return 0; } // not available
 };
 
-// TODO: MQuadrangle8
-// virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
-// virtual char *getStringForBDF(){ return "CQUAD8"; }
-
 class MTetrahedron : public MElement {
  protected:
   MVertex *_v[4];
@@ -407,7 +433,7 @@ class MTetrahedron : public MElement {
 		 _v[trifaces_tetra[num][1]],
 		 _v[trifaces_tetra[num][2]]);
   }
-  virtual int getTypeForMSH(){ return TET1; }
+  virtual int getTypeForMSH(){ return TET_4; }
   virtual int getTypeForUNV(){ return 111; } // solid linear tetrahedron
   virtual char *getStringForPOS(){ return "SS"; }
   virtual char *getStringForBDF(){ return "CTETRA"; }
@@ -464,10 +490,9 @@ class MTetrahedron10 : public MTetrahedron {
   }
   virtual int getNumEdgeVertices(){ return 6; }
   // TODO: edgeRep, faceRep
-  virtual int getTypeForMSH(){ return TET2; }
+  virtual int getTypeForMSH(){ return TET_10; }
   virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
   virtual char *getStringForPOS(){ return "SS2"; }
-  virtual char *getStringForBDF(){ return 0; } // TODO: "CTETRA"
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -512,7 +537,7 @@ class MHexahedron : public MElement {
 		 _v[quadfaces_hexa[num][2]],
 		 _v[quadfaces_hexa[num][3]]);
   }
-  virtual int getTypeForMSH(){ return HEX1; }
+  virtual int getTypeForMSH(){ return HEX_8; }
   virtual int getTypeForUNV(){ return 115; } // solid linear brick
   virtual char *getStringForPOS(){ return "SH"; }
   virtual char *getStringForBDF(){ return "CHEXA"; }
@@ -540,6 +565,56 @@ class MHexahedron : public MElement {
   }
 };
 
+class MHexahedron20 : public MHexahedron {
+ protected:
+  MVertex *_vs[12];
+ public :
+  MHexahedron20(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, 
+		MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
+		MVertex *v10, MVertex *v11, MVertex *v12, MVertex *v13, MVertex *v14,
+		MVertex *v15, MVertex *v16, MVertex *v17, MVertex *v18, MVertex *v19,
+		int num=0, int part=0) 
+    : MHexahedron(v0, v1, v2, v3, v4, v5, v6, v7, num, part)
+  {
+    _vs[0] = v8; _vs[1] = v9; _vs[2] = v10; _vs[3] = v11; _vs[4] = v12; 
+    _vs[5] = v13; _vs[6] = v14; _vs[7] = v15; _vs[8] = v16; _vs[9] = v17; 
+    _vs[10] = v18; _vs[11] = v19; 
+  }
+  MHexahedron20(std::vector<MVertex*> &v, int num=0, int part=0) 
+    : MHexahedron(v, num, part)
+  {
+    for(int i = 0; i < 12; i++) _vs[i] = v[8 + i];
+  }
+  ~MHexahedron20(){}
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 20; }
+  virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+  virtual MVertex *getVertexUNV(int num)
+  {
+    static int map[20] = {0, 8, 1, 11, 2, 13, 3, 9, 10, 12, 
+			  14, 15, 4, 16, 5, 18, 6, 19, 7, 17};
+    return getVertex(map[num]); 
+  }
+  virtual int getNumEdgeVertices(){ return 12; }
+  // TODO: edgeRep, faceRep
+  virtual int getTypeForMSH(){ return HEX_20; }
+  virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
+  virtual void setVolumePositive()
+  {
+    if(getVolumeSign() < 0){
+      MVertex *tmp;
+      tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
+      tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
+      MVertex *old[12];
+      for(int i = 0; i < 12; i++) old[i] = _vs[i];
+      _vs[0] = old[3]; _vs[1] = old[5]; _vs[2] = old[6];
+      _vs[3] = old[0]; _vs[4] = old[4]; _vs[5] = old[1];
+      _vs[6] = old[2]; _vs[7] = old[7]; _vs[8] = old[10];
+      _vs[9] = old[11]; _vs[10] = old[8]; _vs[11] = old[9];
+    }
+  }
+};
+
 class MHexahedron27 : public MHexahedron {
  protected:
   MVertex *_vs[19];
@@ -570,10 +645,8 @@ class MHexahedron27 : public MHexahedron {
   virtual int getNumFaceVertices(){ return 6; }
   virtual int getNumVolumeVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  virtual int getTypeForMSH(){ return HEX2; }
-  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual int getTypeForMSH(){ return HEX_27; }
   virtual char *getStringForPOS(){ return "SH2"; }
-  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -590,10 +663,6 @@ class MHexahedron27 : public MHexahedron {
   }
 };
 
-// TODO: MHexahedron20
-// virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
-// virtual char *getStringForBDF(){ return "CHEXA"; }
-
 class MPrism : public MElement {
  protected:
   MVertex *_v[6];
@@ -632,7 +701,7 @@ class MPrism : public MElement {
 		   _v[quadfaces_prism[num - 2][2]],
 		   _v[quadfaces_prism[num - 2][3]]);
   }
-  virtual int getTypeForMSH(){ return PRI1; }
+  virtual int getTypeForMSH(){ return PRI_6; }
   virtual int getTypeForUNV(){ return 112; } // solid linear wedge
   virtual char *getStringForPOS(){ return "SI"; }
   virtual char *getStringForBDF(){ return "CPENTA"; }
@@ -660,6 +729,50 @@ class MPrism : public MElement {
   }
 };
 
+class MPrism15 : public MPrism {
+ protected:
+  MVertex *_vs[9];
+ public :
+  MPrism15(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, 
+	   MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
+	   MVertex *v10, MVertex *v11, MVertex *v12, MVertex *v13, MVertex *v14,
+	   int num=0, int part=0) 
+    : MPrism(v0, v1, v2, v3, v4, v5, num, part)
+  {
+    _vs[0] = v6; _vs[1] = v7; _vs[2] = v8; _vs[3] = v9; _vs[4] = v10; 
+    _vs[5] = v11; _vs[6] = v12; _vs[7] = v13; _vs[8] = v14;
+  }
+  MPrism15(std::vector<MVertex*> &v, int num=0, int part=0) 
+    : MPrism(v, num, part)
+  {
+    for(int i = 0; i < 9; i++) _vs[i] = v[6 + i];
+  }
+  ~MPrism15(){}
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 15; }
+  virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
+  virtual MVertex *getVertexUNV(int num)
+  {
+    static int map[15] = {0, 6, 1, 9, 2, 7, 8, 10, 11, 3, 12, 4, 14, 5, 13};
+    return getVertex(map[num]); 
+  }
+  virtual int getNumEdgeVertices(){ return 9; }
+  // TODO: edgeRep, faceRep
+  virtual int getTypeForMSH(){ return PRI_15; }
+  virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
+  virtual void setVolumePositive()
+  {
+    if(getVolumeSign() < 0){
+      MVertex *tmp;
+      tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
+      tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
+      tmp = _vs[1]; _vs[1] = _vs[3]; _vs[3] = tmp;
+      tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp;
+      tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp;
+    }
+  }
+};
+
 class MPrism18 : public MPrism {
  protected:
   MVertex *_vs[12];
@@ -686,10 +799,8 @@ class MPrism18 : public MPrism {
   virtual int getNumEdgeVertices(){ return 9; }
   virtual int getNumFaceVertices(){ return 3; }
   // TODO: edgeRep, faceRep
-  virtual int getTypeForMSH(){ return PRI2; }
-  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual int getTypeForMSH(){ return PRI_18; }
   virtual char *getStringForPOS(){ return "SI2"; }
-  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -703,10 +814,6 @@ class MPrism18 : public MPrism {
   }
 };
 
-// TODO: MPrism15
-// virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
-// virtual char *getStringForBDF(){ return "CPENTA"; }
-
 class MPyramid : public MElement {
  protected:
   MVertex *_v[5];
@@ -744,8 +851,7 @@ class MPyramid : public MElement {
 		   _v[quadfaces_pyramid[num - 4][2]],
 		   _v[quadfaces_pyramid[num - 4][3]]);
   }
-  virtual int getTypeForMSH(){ return PYR1; }
-  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual int getTypeForMSH(){ return PYR_5; }
   virtual char *getStringForPOS(){ return "SY"; }
   virtual char *getStringForBDF(){ return "CPYRAM"; }
   virtual int getVolumeSign()
@@ -771,6 +877,42 @@ class MPyramid : public MElement {
   }
 };
 
+class MPyramid13 : public MPyramid {
+ protected:
+  MVertex *_vs[8];
+ public :
+  MPyramid13(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, 
+	     MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, MVertex *v9,
+	     MVertex *v10, MVertex *v11, MVertex *v12, int num=0, int part=0) 
+    : MPyramid(v0, v1, v2, v3, v4, num, part)
+  {
+    _vs[0] = v5; _vs[1] = v6; _vs[2] = v7; _vs[3] = v8; _vs[4] = v9; 
+    _vs[5] = v10; _vs[6] = v11; _vs[7] = v12;
+  }
+  MPyramid13(std::vector<MVertex*> &v, int num=0, int part=0) 
+    : MPyramid(v, num, part)
+  {
+    for(int i = 0; i < 8; i++) _vs[i] = v[5 + i];
+  }
+  ~MPyramid13(){}
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 13; }
+  virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
+  virtual int getNumEdgeVertices(){ return 8; }
+  // TODO: edgeRep, faceRep
+  virtual int getTypeForMSH(){ return PYR_13; }
+  virtual void setVolumePositive()
+  {
+    if(getVolumeSign() < 0){
+      MVertex *tmp;
+      tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
+      tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
+      tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp;
+      tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp;
+    }
+  }
+};
+
 class MPyramid14 : public MPyramid {
  protected:
   MVertex *_vs[9];
@@ -796,10 +938,8 @@ class MPyramid14 : public MPyramid {
   virtual int getNumEdgeVertices(){ return 8; }
   virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  virtual int getTypeForMSH(){ return PYR2; }
-  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual int getTypeForMSH(){ return PYR_14; }
   virtual char *getStringForPOS(){ return "SY2"; }
-  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index c1d3823cf33fca4cb46bfab0469987d028dea885..3a5579763e4f475787dc7574a057a660ffbd5621 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.96 2006-09-03 07:44:10 geuzaine Exp $
+// $Id: Generator.cpp,v 1.97 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -101,9 +101,6 @@ void GetStatistics(double stat[50], double quality[3][100])
   stat[14] = CTX.mesh_timer[1];
   stat[15] = CTX.mesh_timer[2];
   
-  // FIXME:
-  //stat[16] = numOrder2Vertices; 
-  
   if(quality){
     for(int i = 0; i < 3; i++)
       for(int j = 0; j < 100; j++)
@@ -436,7 +433,7 @@ void mai3d(int ask)
 
   // Create second order elements
   if(GMODEL->getMeshStatus() && CTX.mesh.order == 2) 
-    Degre2();
+    Degre2(CTX.mesh.second_order_linear, CTX.mesh.second_order_incomplete);
 
   // Partition
   if(GMODEL->getMeshStatus() > 1 && CTX.mesh.nbPartitions != 1)
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 8ae5cf31210a389b0cb3bb8c38efceb7e05b58c9..6a42091f517a6ee45a767fb8d594639e4650f016 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -397,6 +397,6 @@ int Recombine_All(Mesh *M);
 void ApplyLcFactor();
 
 void Degre1();
-void Degre2();
+void Degre2(bool linear=true, bool incomplete=false);
 
 #endif
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index 08505757027c1aaa5cdfa0290c99661f4a355048..54a3b93fb9828a3c0520d6642eb9cf9128131a43 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.45 2006-09-07 19:45:15 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.46 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -206,7 +206,7 @@ void setSecondOrder(GEdge *ge,
 void setSecondOrder(GFace *gf,
 		    std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
 		    std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
-		    bool linear)
+		    bool linear, bool incomplete)
 {
   std::vector<MTriangle*> triangles2;
   for(unsigned int i = 0; i < gf->triangles.size(); i++){
@@ -225,10 +225,17 @@ void setSecondOrder(GFace *gf,
     MQuadrangle *q = gf->quadrangles[i];
     std::vector<MVertex*> ve, vf;
     getEdgeVertices(gf, q, ve, edgeVertices, linear);
-    getFaceVertices(gf, q, vf, faceVertices, linear);
-    quadrangles2.push_back
-      (new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),
-			q->getVertex(3), ve[0], ve[1], ve[2], ve[3], vf[0]));
+    if(incomplete){
+      quadrangles2.push_back
+	(new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2),
+			  q->getVertex(3), ve[0], ve[1], ve[2], ve[3]));
+    }
+    else{
+      getFaceVertices(gf, q, vf, faceVertices, linear);
+      quadrangles2.push_back
+	(new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),
+			  q->getVertex(3), ve[0], ve[1], ve[2], ve[3], vf[0]));
+    }
     delete q;
   }
   gf->quadrangles = quadrangles2;
@@ -239,7 +246,7 @@ void setSecondOrder(GFace *gf,
 void setSecondOrder(GRegion *gr,
 		    std::map<std::pair<MVertex*,MVertex*>, MVertex* > &edgeVertices,
 		    std::map<std::vector<MVertex*>, MVertex* > &faceVertices,
-		    bool linear)
+		    bool linear, bool incomplete)
 {
   std::vector<MTetrahedron*> tetrahedra2;
   for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
@@ -258,16 +265,26 @@ void setSecondOrder(GRegion *gr,
     MHexahedron *h = gr->hexahedra[i];
     std::vector<MVertex*> ve, vf;
     getEdgeVertices(gr, h, ve, edgeVertices, linear);
-    getFaceVertices(gr, h, vf, faceVertices, linear);
-    SPoint3 pc = h->barycenter();
-    MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
-    gr->mesh_vertices.push_back(v);
-    hexahedra2.push_back
-      (new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2), 
-			 h->getVertex(3), h->getVertex(4), h->getVertex(5), 
-			 h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2], 
-			 ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10], 
-			 ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v));
+    if(incomplete){
+      hexahedra2.push_back
+	(new MHexahedron20(h->getVertex(0), h->getVertex(1), h->getVertex(2), 
+			   h->getVertex(3), h->getVertex(4), h->getVertex(5), 
+			   h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2], 
+			   ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10], 
+			   ve[11]));
+    }
+    else{
+      getFaceVertices(gr, h, vf, faceVertices, linear);
+      SPoint3 pc = h->barycenter();
+      MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
+      gr->mesh_vertices.push_back(v);
+      hexahedra2.push_back
+	(new MHexahedron27(h->getVertex(0), h->getVertex(1), h->getVertex(2), 
+			   h->getVertex(3), h->getVertex(4), h->getVertex(5), 
+			   h->getVertex(6), h->getVertex(7), ve[0], ve[1], ve[2], 
+			   ve[3], ve[4], ve[5], ve[6], ve[7], ve[8], ve[9], ve[10], 
+			   ve[11], vf[0], vf[1], vf[2], vf[3], vf[4], vf[5], v));
+    }
     delete h;
   }
   gr->hexahedra = hexahedra2;
@@ -277,12 +294,20 @@ void setSecondOrder(GRegion *gr,
     MPrism *p = gr->prisms[i];
     std::vector<MVertex*> ve, vf;
     getEdgeVertices(gr, p, ve, edgeVertices, linear);
-    getFaceVertices(gr, p, vf, faceVertices, linear);
-    prisms2.push_back
-      (new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
-		    p->getVertex(3), p->getVertex(4), p->getVertex(5), 
-		    ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7],
-		    ve[8], vf[0], vf[1], vf[2]));
+    if(incomplete){
+      prisms2.push_back
+	(new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+		      p->getVertex(3), p->getVertex(4), p->getVertex(5), 
+		      ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8]));
+    }
+    else{
+      getFaceVertices(gr, p, vf, faceVertices, linear);
+      prisms2.push_back
+	(new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+		      p->getVertex(3), p->getVertex(4), p->getVertex(5), 
+		      ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
+		      vf[0], vf[1], vf[2]));
+    }
     delete p;
   }
   gr->prisms = prisms2;
@@ -292,11 +317,19 @@ void setSecondOrder(GRegion *gr,
     MPyramid *p = gr->pyramids[i];
     std::vector<MVertex*> ve, vf;
     getEdgeVertices(gr, p, ve, edgeVertices, linear);
-    getFaceVertices(gr, p, vf, faceVertices, linear);
-    pyramids2.push_back
-      (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
-		      p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
-		      ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
+    if(incomplete){
+      pyramids2.push_back
+	(new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+			p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
+			ve[3], ve[4], ve[5], ve[6], ve[7]));
+    }
+    else{
+      getFaceVertices(gr, p, vf, faceVertices, linear);
+      pyramids2.push_back
+	(new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+			p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
+			ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
+    }
     delete p;
   }
   gr->pyramids = pyramids2;
@@ -364,30 +397,34 @@ void Degre1()
     removeSecondOrderVertices(*it);
 }
 
-void Degre2()
+void Degre2(bool linear, bool incomplete)
 {
   Msg(STATUS1, "Meshing second order...");
   double t1 = Cpu();
 
-  bool linear = true;
-  //bool linear = false;
+  // This routine replaces all the elements in the mesh with second
+  // order elements by creating unique nodes on the edges/faces of the
+  // mesh:
+  // - If linear is set to true, new vertices are created by linear
+  //   interpolation between existing ones. If not, new vertices are
+  //   created on the exact geometry, provided that the geometrical
+  //   edges/faces are discretized with 1D/2D elements. (I.e., if
+  //   there are only 3D elements in the mesh then any new nodes on
+  //   the boundary will always be created by linear interpolation,
+  //   whether linear is set or not.)
+  // - If incomplete is set to true, we only create new vertices on 
+  //   edges (creating 8-node quads, 20-node hexas, etc., instead of
+  //   9-node quads, 27-node hexas, etc.)
 
   std::map<std::pair<MVertex*,MVertex*>, MVertex* > edgeVertices;
   std::map<std::vector<MVertex*>, MVertex* > faceVertices;
 
-  // replace all elements with second order elements by creating
-  // unique nodes on the edges/faces of the mesh. (To generate nodes
-  // on the exact geometrical edges/faces this assumes that the
-  // geometrical edges/faces are discretized with 1D/2D
-  // elements. I.e., if there are only 3D elements in the mesh then
-  // any new nodes on the boundary will simply be added by linear
-  // interpolation.)
   for(GModel::eiter it = GMODEL->firstEdge(); it != GMODEL->lastEdge(); ++it)
     setSecondOrder(*it, edgeVertices, linear);
   for(GModel::fiter it = GMODEL->firstFace(); it != GMODEL->lastFace(); ++it)
-    setSecondOrder(*it, edgeVertices, faceVertices, linear);
+    setSecondOrder(*it, edgeVertices, faceVertices, linear, incomplete);
   for(GModel::riter it = GMODEL->firstRegion(); it != GMODEL->lastRegion(); ++it)
-    setSecondOrder(*it, edgeVertices, faceVertices, linear);
+    setSecondOrder(*it, edgeVertices, faceVertices, linear, incomplete);
 
   double t2 = Cpu();
   Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4f204b1f336f0bc50a0fe9c5490fb69c460e02e4..dcbc2c7585363f423802fe8c045bcda67bc7a414 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.11 2006-09-05 21:37:59 remacle Exp $
+// $Id: meshGFace.cpp,v 1.12 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -294,13 +294,10 @@ void OptimizeMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 
 void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 {
-  int NUMP = 0;
   int IT =0;
 
   int MAXNP = m.MAXPOINTNUMBER;
 
-  clock_t t1 = clock();
-
   // computecharacteristic lengths using mesh edge spacing
   
   std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
@@ -359,7 +356,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
 	  ++it;
 	}
 
-      //      Msg(STATUS2," %d triangles : conv %g -> %g (%g sec)",m.triangles.size(),maxL,1.5,(double)(clock()-t1)/CLOCKS_PER_SEC);
+
       if ((minL > 0.4 && maxL < 1.5) || IT > NIT)break;
 
 
@@ -626,7 +623,7 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
   if (e)e->g = g;
   else throw;
 
-  for (int i=1;i<ge->mesh_vertices.size();i++)
+  for (unsigned int i=1;i<ge->mesh_vertices.size();i++)
     {
       vstart = ge->mesh_vertices[i-1];
       vend   = ge->mesh_vertices[i];
@@ -796,7 +793,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
   for(int i = 0; i < doc.numPoints; i++) 
     {
       MVertex *here = (MVertex*)doc.points[i].data;      
-      BDS_Point *bds_p = m->add_point ( here->getNum(), here->x(), here->y(), gf);
+      m->add_point ( here->getNum(), here->x(), here->y(), gf);
     }
   for(int i = 0; i < doc.numTriangles; i++) 
     {
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index 1a191019ebd5f8d733b29ed36c590611babd526d..34c08972d291e0eb5d881cdb618b93898b4b7ece 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.3 2006-09-07 16:03:32 remacle Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.4 2006-09-08 02:39:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -40,7 +40,7 @@ int MeshTransfiniteSurface( GFace *gf)
   std::vector <MVertex *> d_vertices;
   std::vector <int> indices;
 
-  for (int i=0;i<gf->meshAttributes.corners.size();i++)
+  for (unsigned int i=0;i<gf->meshAttributes.corners.size();i++)
     corners.push_back(gf->meshAttributes.corners[i]->mesh_vertices[0]);
 
   computeEdgeLoops (gf,d_vertices, indices);
@@ -51,10 +51,10 @@ int MeshTransfiniteSurface( GFace *gf)
     Msg (GERROR,"Surface %d is transfinite but has %d holes",gf->tag(),indices.size()-2);
 
   std::vector <MVertex *> m_vertices;
-  int I;
+  unsigned int I;
   for (I=0;I<d_vertices.size();I++)
     if(d_vertices[I] == corners[0])break;
-  for (int j=0;j<d_vertices.size();j++)
+  for (unsigned int j=0;j<d_vertices.size();j++)
     m_vertices.push_back(d_vertices[(I+j)%d_vertices.size()]);
   
   int iCorner = 0;
@@ -64,7 +64,7 @@ int MeshTransfiniteSurface( GFace *gf)
 
   std::map<std::pair<int,int> , MVertex*> tab;
 
-  for (int i=0;i<m_vertices.size();i++)
+  for (unsigned int i=0;i<m_vertices.size();i++)
     {
       MVertex *v = m_vertices[i];
       if (v==corners[0]||v==corners[1]||v==corners[2]|| (corners.size()==4 && v==corners[3]))
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 1f7fe659d96c285bdf6040ec412e7ede58c75cbf..92d678ffe16a262866f3eb4f1203a3235af5e1b2 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -16,7 +16,7 @@ void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex
   while (it != faces.end())
     {
       GFace *gf = (*it);      
-      for (int i = 0; i< gf->triangles.size(); i++)
+      for (unsigned int i = 0; i< gf->triangles.size(); i++)
 	{
 	  MTriangle *t = gf->triangles[i];
 	  allBoundingVertices.insert (t->getVertex(0));
@@ -77,7 +77,7 @@ void buildTetgenStructure (  GRegion *gr, tetgenio &in, std::vector<MVertex*> &
   while (it != faces.end())
     {
       GFace *gf = (*it);
-      for (int i = 0; i< gf->triangles.size(); i++)
+      for (unsigned int i = 0; i< gf->triangles.size(); i++)
 	{
 	  MTriangle *t = gf->triangles[i];
 	  tetgenio::facet *f = &in.facetlist[I];
@@ -101,7 +101,6 @@ void buildTetgenStructure (  GRegion *gr, tetgenio &in, std::vector<MVertex*> &
 
 void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MVertex*> & numberedV)
 {
-  int I = numberedV.size() + 1;
   for (int i = 0; i < out.numberofpoints; i++) 
     {
       MVertex *v = new MVertex (out.pointlist[i * 3 + 0],out.pointlist[i * 3 + 1],out.pointlist[i * 3 + 2]); 
@@ -159,7 +158,7 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
   while (it != faces.end())
     {
       GFace *gf = (*it);      
-      for (int i = 0; i< gf->triangles.size(); i++)
+      for (unsigned int i = 0; i< gf->triangles.size(); i++)
 	{
 	  MTriangle *t = gf->triangles[i];
 	  int tmp[3];
@@ -173,7 +172,7 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
 
   if (importVolumeMesh)
     {
-      for (int i = 0; i< gr->tetrahedra.size(); i++)
+      for (unsigned int i = 0; i< gr->tetrahedra.size(); i++)
 	{
 	  MTetrahedron *t = gr->tetrahedra[i];
 	  int tmp[4];
@@ -298,7 +297,7 @@ void meshNormalsPointOutOfTheRegion (GRegion *gr)
     {
       GFace *gf = (*it);      
       int nb_intersect = 0;
-      for (int i = 0; i< gf->triangles.size(); i++)
+      for (unsigned int i = 0; i< gf->triangles.size(); i++)
 	{
 	  MTriangle *t = gf->triangles[i];
 	  double X[3] = {t->getVertex(0)->x(),t->getVertex(1)->x(),t->getVertex(2)->x()};
@@ -320,7 +319,7 @@ void meshNormalsPointOutOfTheRegion (GRegion *gr)
 	  while (it_b != faces.end())
 	    {
 	      GFace *gf_b = (*it_b);
-	      for (int i_b = 0; i_b< gf_b->triangles.size(); i_b++)
+	      for (unsigned int i_b = 0; i_b< gf_b->triangles.size(); i_b++)
 		{
 		  MTriangle *t_b = gf_b->triangles[i_b];
 		  if (t_b != t)
@@ -341,7 +340,7 @@ void meshNormalsPointOutOfTheRegion (GRegion *gr)
 
       if (nb_intersect % 2 == 1) // odd nb of intersections: the normal points inside the region 
 	{
-	  for (int i = 0; i< gf->triangles.size(); i++)
+	  for (unsigned int i = 0; i< gf->triangles.size(); i++)
 	    {
 	      gf->triangles[i]->revert();
 	    }
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 67c21439f8b5095b6e5be577eff33c2a508a2d1c..d69b0c796327dfb0c88ee21719e8e8b690157b77 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,11 +1,12 @@
-$Id: VERSIONS,v 1.365 2006-09-03 07:44:11 geuzaine Exp $
+$Id: VERSIONS,v 1.366 2006-09-08 02:39:43 geuzaine Exp $
 
 2.0: new geometry and mesh databases; complete rewrite of geometry and
 mesh drawing code; complete rewrite of the input/output code (with new
 native binary MSH format and full support for import/export of I-deas
-UNV, Nastran BDF, STL, Medit MESH and VRML 1.0 files); new 2D mesh
-algorithm; option changes in the interface are now applied
-instantenously; lots of improvements all over the place.
+UNV, Nastran BDF, STL, Medit MESH and VRML 1.0 files); added support
+for incomplete second order elements; new 2D mesh algorithm; option
+changes in the interface are now applied instantenously; lots of
+improvements all over the place.
 
 1.66: added support for offscreen rendering using OSMesa; added
 support for SVG output;
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 4e10f5293bcd6a41e36ff61c6f77d9a05da9199f..1152ac1d12707343d99e9667d22791a2604f1473 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.215 2006-09-07 02:56:53 geuzaine Exp $
+@c $Id: gmsh.texi,v 1.216 2006-09-08 02:39:43 geuzaine Exp $
 @c
 @c Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 @c
@@ -2929,42 +2929,54 @@ ordered) way.
 defines the geometrical type of the @var{n}-th element:
 @table @code
 @item 1
-Line (2 nodes).
+2-node line.
 @item 2
-Triangle (3 nodes).
+3-node triangle.
 @item 3
-Quadrangle (4 nodes).
+4-node quadrangle.
 @item 4
-Tetrahedron (4 nodes).
+4-node tetrahedron.
 @item 5
-Hexahedron (8 nodes).
+8-node hexahedron.
 @item 6
-Prism (6 nodes).
+6-node prism.
 @item 7
-Pyramid (5 nodes).
+5-node pyramid.
 @item 8
-Second order line (3 nodes: 2 associated with the vertices and 1 with the
-edge).
+3-node second order line (2 nodes associated with the vertices and 1
+with the edge).
 @item 9
-Second order triangle (6 nodes: 3 associated with the vertices and 3 with
-the edges).
+6-node second order triangle (3 nodes associated with the vertices and 3
+with the edges).
 @item 10
-Second order quadrangle (9 nodes: 4 associated with the vertices, 4 with the
-edges and 1 with the face).
+9-node second order quadrangle (4 nodes associated with the vertices, 4
+with the edges and 1 with the face).
 @item 11
-Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
-with the edges).
+10-node second order tetrahedron (4 nodes associated with the vertices
+and 6 with the edges).
 @item 12
-Second order hexahedron (27 nodes: 8 associated with the vertices, 12 with
-the edges, 6 with the faces and 1 with the volume).
+27-node second order hexahedron (8 nodes associated with the vertices,
+12 with the edges, 6 with the faces and 1 with the volume).
 @item 13
-Second order prism (18 nodes: 6 associated with the vertices, 9 with the
-edges and 3 with the quadrangular faces).
+18-node second order prism (6 nodes associated with the vertices, 9 with
+the edges and 3 with the quadrangular faces).
 @item 14
-Second order pyramid (14 nodes: 5 associated with the vertices, 8 with the
-edges and 1 with the quadrangular face).
+14-node second order pyramid (5 nodes associated with the vertices, 8
+with the edges and 1 with the quadrangular face).
 @item 15
-Point (1 node).
+1-node point.
+@item 16
+8-node second order quadrangle (4 nodes associated with the vertices and
+4 with the edges).
+@item 17
+20-node second order hexahedron (8 nodes associated with the vertices
+and 12 with the edges).
+@item 18
+15-node second order prism (6 nodes associated with the vertices and 9
+with the edges).
+@item 19
+13-node second order pyramid (5 nodes associated with the vertices and 8
+with the edges).
 @end table
 See below for the ordering of the nodes.
 
@@ -3056,42 +3068,54 @@ ordered) way.
 defines the geometrical type of the @var{n}-th element:
 @table @code
 @item 1
-Line (2 nodes).
+2-node line.
 @item 2
-Triangle (3 nodes).
+3-node triangle.
 @item 3
-Quadrangle (4 nodes).
+4-node quadrangle.
 @item 4
-Tetrahedron (4 nodes).
+4-node tetrahedron.
 @item 5
-Hexahedron (8 nodes).
+8-node hexahedron.
 @item 6
-Prism (6 nodes).
+6-node prism.
 @item 7
-Pyramid (5 nodes).
+5-node pyramid.
 @item 8
-Second order line (3 nodes: 2 associated with the vertices and 1 with the
-edge).
+3-node second order line (2 nodes associated with the vertices and 1
+with the edge).
 @item 9
-Second order triangle (6 nodes: 3 associated with the vertices and 3 with
-the edges).
+6-node second order triangle (3 nodes associated with the vertices and 3
+with the edges).
 @item 10
-Second order quadrangle (9 nodes: 4 associated with the vertices, 4 with the
-edges and 1 with the face).
+9-node second order quadrangle (4 nodes associated with the vertices, 4
+with the edges and 1 with the face).
 @item 11
-Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
-with the edges).
+10-node second order tetrahedron (4 nodes associated with the vertices
+and 6 with the edges).
 @item 12
-Second order hexahedron (27 nodes: 8 associated with the vertices, 12 with
-the edges, 6 with the faces and 1 with the volume).
+27-node second order hexahedron (8 nodes associated with the vertices,
+12 with the edges, 6 with the faces and 1 with the volume).
 @item 13
-Second order prism (18 nodes: 6 associated with the vertices, 9 with the
-edges and 3 with the quadrangular faces).
+18-node second order prism (6 nodes associated with the vertices, 9 with
+the edges and 3 with the quadrangular faces).
 @item 14
-Second order pyramid (14 nodes: 5 associated with the vertices, 8 with the
-edges and 1 with the quadrangular face).
+14-node second order pyramid (5 nodes associated with the vertices, 8
+with the edges and 1 with the quadrangular face).
 @item 15
-Point (1 node).
+1-node point.
+@item 16
+8-node second order quadrangle (4 nodes associated with the vertices and
+4 with the edges).
+@item 17
+20-node second order hexahedron (8 nodes associated with the vertices
+and 12 with the edges).
+@item 18
+15-node second order prism (6 nodes associated with the vertices and 9
+with the edges).
+@item 19
+13-node second order pyramid (5 nodes associated with the vertices and 8
+with the edges).
 @end table
 See below for the ordering of the nodes.
 
diff --git a/doc/texinfo/opt_general.texi b/doc/texinfo/opt_general.texi
index d9e39e8ce25823aabedd6022948d71ff8c139f82..c77f751a372584cd3b4ad8aefc93d02c88cb23ca 100644
--- a/doc/texinfo/opt_general.texi
+++ b/doc/texinfo/opt_general.texi
@@ -866,7 +866,7 @@ Saved in: @code{General.OptionsFileName}
 
 @item General.Color.Foreground
 Foreground color@*
-Default value: @code{@{255,255,255@}}@*
+Default value: @code{@{170,170,170@}}@*
 Saved in: @code{General.OptionsFileName}
 
 @item General.Color.Text
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index 945c9d3ee4176ec5bd03d7ace6b894c7cc91686e..14eb045957472cf07edd8e4297be9b094a6a9fcf 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -34,6 +34,11 @@ Maximum ratio of two consecutive edge lengths@*
 Default value: @code{0.9}@*
 Saved in: @code{General.OptionsFileName}
 
+@item Mesh.BdfFieldFormat
+Field format for Nastran BDF files (0=free, 1=small, 2=large)@*
+Default value: @code{1}@*
+Saved in: @code{General.OptionsFileName}
+
 @item Mesh.CharacteristicLengthFactor
 Factor applied to all characteristic lengths (and background meshes)@*
 Default value: @code{1}@*
@@ -115,7 +120,7 @@ Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
 @item Mesh.Format
-Mesh output format (1=msh, 2=unv, 3=gref, 19=vrml)@*
+Mesh output format (1=msh, 2=unv, 19=vrml, 27=stl, 30=mesh, 31=bdf)@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
@@ -186,7 +191,7 @@ Saved in: @code{General.OptionsFileName}
 
 @item Mesh.MshFileVersion
 Version of the MSH file format to use@*
-Default value: @code{1}@*
+Default value: @code{2}@*
 Saved in: @code{General.OptionsFileName}
 
 @item Mesh.NbElemsPerRadiusOfCurv
@@ -239,6 +244,11 @@ Optimize the mesh using Netgen to improve the quality of tetrahedral elements@*
 Default value: @code{0}@*
 Saved in: @code{General.OptionsFileName}
 
+@item Mesh.Partitioning
+Number of partitions applied to the final mesh@*
+Default value: @code{1}@*
+Saved in: @code{General.OptionsFileName}
+
 @item Mesh.Points
 Display mesh vertices (nodes)?@*
 Default value: @code{0}@*
@@ -334,13 +344,18 @@ Global scaling factor applied to the saved mesh@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
-@item Mesh.Smoothing
-Number of smoothing steps applied to the final mesh@*
+@item Mesh.SecondOrderIncomplete
+Create incomplete second order elements? (8-node quads, 20-node hexas, etc.)@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
-@item Mesh.Partitioning
-Number of partitions applied to the final mesh@*
+@item Mesh.SecondOrderLinear
+Should second order vertices simply be created by linear interpolation?@*
+Default value: @code{1}@*
+Saved in: @code{General.OptionsFileName}
+
+@item Mesh.Smoothing
+Number of smoothing steps applied to the final mesh@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}