diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index 3765c08f9f5ac493cb2d6a07f5a3a8d55a2171bc..7d3416e3ca315e665dacbb8cb10e0f1db6a004c4 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.47 2006-09-08 13:54:27 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.48 2006-09-08 17:45:51 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -1596,44 +1596,45 @@ int GModel::readBDF(const std::string &name)
     if(!fgets(buffer, sizeof(buffer), fp)) break;
     if(buffer[0] != '$'){ // skip comments
       int num, region;
-      std::vector<MVertex*> n;
+      std::vector<MVertex*> vertices;
       if(!strncmp(buffer, "CBAR", 4)){
-	if(readElementBDF(fp, buffer, 4, 2, &num, &region, n, vertexMap))
-	  elements[0][region].push_back(new MLine(n, num));
+	if(readElementBDF(fp, buffer, 4, 2, &num, &region, vertices, vertexMap))
+	  elements[0][region].push_back(new MLine(vertices, num));
       }
       else if(!strncmp(buffer, "CTRIA3", 6)){
-	if(readElementBDF(fp, buffer, 6, 3, &num, &region, n, vertexMap))
-	  elements[1][region].push_back(new MTriangle(n, num));
+	if(readElementBDF(fp, buffer, 6, 3, &num, &region, vertices, vertexMap))
+	  elements[1][region].push_back(new MTriangle(vertices, num));
       }
       else if(!strncmp(buffer, "CTRIA6", 6)){
-	// not sure about the node ordering
-	if(readElementBDF(fp, buffer, 6, 6, &num, &region, n, vertexMap))
-	  elements[1][region].push_back(new MTriangle6(n, num));
+	if(readElementBDF(fp, buffer, 6, 6, &num, &region, vertices, vertexMap))
+	  elements[1][region].push_back(new MTriangle6(vertices, num));
       }
       else if(!strncmp(buffer, "CQUAD4", 6)){
-	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertexMap))
-	  elements[2][region].push_back(new MQuadrangle(n, num));
+	if(readElementBDF(fp, buffer, 6, 4, &num, &region, vertices, vertexMap))
+	  elements[2][region].push_back(new MQuadrangle(vertices, num));
       }
       else if(!strncmp(buffer, "CQUAD8", 6)){
-	// not sure about the node ordering
-	if(readElementBDF(fp, buffer, 6, 8, &num, &region, n, vertexMap))
-	  elements[2][region].push_back(new MQuadrangle8(n, num));
+	if(readElementBDF(fp, buffer, 6, 8, &num, &region, vertices, vertexMap))
+	  elements[2][region].push_back(new MQuadrangle8(vertices, num));
       }
       else if(!strncmp(buffer, "CTETRA", 6)){
-	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertexMap))
-	  elements[3][region].push_back(new MTetrahedron(n, num));
+	// should be generalized to also read 10-node tets
+	if(readElementBDF(fp, buffer, 6, 4, &num, &region, vertices, vertexMap))
+	  elements[3][region].push_back(new MTetrahedron(vertices, num));
       }
       else if(!strncmp(buffer, "CHEXA", 5)){
-	if(readElementBDF(fp, buffer, 5, 8, &num, &region, n, vertexMap))
-	  elements[4][region].push_back(new MHexahedron(n, num));
+	// should be generalized to also read 20-node hexas
+	if(readElementBDF(fp, buffer, 5, 8, &num, &region, vertices, vertexMap))
+	  elements[4][region].push_back(new MHexahedron(vertices, num));
       }
       else if(!strncmp(buffer, "CPENTA", 6)){
-	if(readElementBDF(fp, buffer, 6, 6, &num, &region, n, vertexMap))
-	  elements[5][region].push_back(new MPrism(n, num));
+	// should be generalized to also read 15-node prisms
+	if(readElementBDF(fp, buffer, 6, 6, &num, &region, vertices, vertexMap))
+	  elements[5][region].push_back(new MPrism(vertices, num));
       }
       else if(!strncmp(buffer, "CPYRAM", 6)){
-	if(readElementBDF(fp, buffer, 6, 5, &num, &region, n, vertexMap))
-	  elements[6][region].push_back(new MPyramid(n, num));
+	if(readElementBDF(fp, buffer, 6, 5, &num, &region, vertices, vertexMap))
+	  elements[6][region].push_back(new MPyramid(vertices, num));
       }
     }
   }
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9ab1d5b82075f851333cc88fb05cc6905041ef2f..30d39d882004cb6518b3513c5273a57a9876143d 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.17 2006-09-08 02:39:43 geuzaine Exp $
+// $Id: MElement.cpp,v 1.18 2006-09-08 17:45:51 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -110,12 +110,13 @@ SPoint3 MElement::barycenter()
 void MElement::writeMSH(FILE *fp, double version, bool binary, int num, 
 			int elementary, int physical)
 {
+  int type = getTypeForMSH();
+  if(!type) return;
+
   // if necessary, change the ordering of the vertices to get positive
   // volume
   setVolumePositive();
-
   int n = getNumVertices();
-  int type = getTypeForMSH();
 
   if(!binary){
     fprintf(fp, "%d %d", num ? num : _num, type);
@@ -164,34 +165,34 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
 void MElement::writePOS(FILE *fp, double scalingFactor, int elementary)
 {
   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");
+  if(!str) return;
+
+  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");
 }
 
 void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
@@ -255,24 +256,24 @@ void MElement::writeVRML(FILE *fp)
 void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
 {
   int type = getTypeForUNV();
-  if(type){
-    setVolumePositive();
-    int n = getNumVertices();
-    int physical_property = elementary;
-    int material_property = physical;
-    int color = 7;
-    fprintf(fp, "%10d%10d%10d%10d%10d%10d\n",
-	    num ? num : _num, type, physical_property, material_property, color, n);
-    if(type == 21 || type == 24) // linear beam or parabolic beam
-      fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);
-    for(int k = 0; k < n; k++) {
-      fprintf(fp, "%10d", getVertexUNV(k)->getNum());
-      if(k % 8 == 7)
-	fprintf(fp, "\n");
-    }
-    if(n - 1 % 8 != 7)
+  if(!type) return;
+
+  setVolumePositive();
+  int n = getNumVertices();
+  int physical_property = elementary;
+  int material_property = physical;
+  int color = 7;
+  fprintf(fp, "%10d%10d%10d%10d%10d%10d\n",
+	  num ? num : _num, type, physical_property, material_property, color, n);
+  if(type == 21 || type == 24) // linear beam or parabolic beam
+    fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);
+  for(int k = 0; k < n; k++) {
+    fprintf(fp, "%10d", getVertexUNV(k)->getNum());
+    if(k % 8 == 7)
       fprintf(fp, "\n");
   }
+  if(n - 1 % 8 != 7)
+    fprintf(fp, "\n");
 }
 
 void MElement::writeMESH(FILE *fp, int elementary)
@@ -285,29 +286,37 @@ void MElement::writeMESH(FILE *fp, int elementary)
 void MElement::writeBDF(FILE *fp, int format, int elementary)
 {
   char *str = getStringForBDF();
-  if(str){
-    int n = getNumVertices();
-    if(format == 0){ // free field format
-      fprintf(fp, "%s,%d,%d", str, _num, elementary);
-      for(int i = 0; i < n; i++){
-	if(i != n - 1 && !((i + 3) % 9))
-	  fprintf(fp, ",+E%d\n+E%d", _num, _num);
-	fprintf(fp, ",%d", getVertex(i)->getNum());
+  if(!str) return;
+
+  setVolumePositive();
+  int n = getNumVertices();
+  const char *cont[4] = {"E", "F", "G", "H"};
+  int ncont = 0;
+  
+  if(format == 0){ // free field format
+    fprintf(fp, "%s,%d,%d", str, _num, elementary);
+    for(int i = 0; i < n; i++){
+      fprintf(fp, ",%d", getVertex(i)->getNum());
+      if(i != n - 1 && !((i + 3) % 8)){
+	fprintf(fp, ",+%s%d\n+%s%d", cont[ncont], _num, cont[ncont], _num);
+	ncont++;
       }
-      if(n == 2) // CBAR
-	fprintf(fp, ",0.,0.,0.");
-      fprintf(fp, "\n");
     }
-    else{ // small or large field format
-      fprintf(fp, "%-8s%-8d%-8d", str, _num, elementary);
-      for(int i = 0; i < n; i++){
-	if(i != n - 1 && !((i + 3) % 9))
-	  fprintf(fp, "+E%-6d\n+E%-6d", _num, _num);
-	fprintf(fp, "%-8d", getVertex(i)->getNum());
+    if(n == 2) // CBAR
+      fprintf(fp, ",0.,0.,0.");
+    fprintf(fp, "\n");
+  }
+  else{ // small or large field format
+    fprintf(fp, "%-8s%-8d%-8d", str, _num, elementary);
+    for(int i = 0; i < n; i++){
+      fprintf(fp, "%-8d", getVertex(i)->getNum());
+      if(i != n - 1 && !((i + 3) % 8)){
+	fprintf(fp, "+%s%-6d\n+%s%-6d", cont[ncont], _num, cont[ncont], _num);
+	ncont++;
       }
-      if(n == 2) // CBAR
-	fprintf(fp, "%-8s%-8s%-8s", "0.", "0.", "0.");
-      fprintf(fp, "\n");
     }
+    if(n == 2) // CBAR
+      fprintf(fp, "%-8s%-8s%-8s", "0.", "0.", "0.");
+    fprintf(fp, "\n");
   }
 }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index d094b4e9fea73636171c69f2cc4b229a98a571b7..3e3886a2ee6ea3740b5a2eeb7627923b5407b436 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -92,6 +92,9 @@ class MElement
   // get the vertex using the I-deas UNV ordering
   virtual MVertex *getVertexUNV(int num){ return getVertex(num); }
 
+  // get the vertex using the Nastran BDF ordering
+  virtual MVertex *getVertexBDF(int num){ return getVertex(num); }
+
   // get the number of vertices associated with edges, faces and
   // volumes (nonzero only for higher order elements)
   virtual int getNumEdgeVertices(){ return 0; }
@@ -142,11 +145,11 @@ class MElement
   virtual void writeMESH(FILE *fp, int elementary=1);
   virtual void writeBDF(FILE *fp, int format=0, int elementary=1);
 
-  // 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; }
+  // info for specific IO formats
+  virtual int getTypeForMSH() = 0;
+  virtual int getTypeForUNV() = 0;
+  virtual char *getStringForPOS() = 0;
+  virtual char *getStringForBDF() = 0;
 };
 
 class MLine : public MElement {
@@ -214,6 +217,7 @@ class MLine3 : public MLine {
   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 {
@@ -374,6 +378,7 @@ class MQuadrangle8 : public MQuadrangle {
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return QUA_8; }
   virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
+  virtual char *getStringForPOS(){ return 0; } // not available
   virtual char *getStringForBDF(){ return "CQUAD8"; }
 };
 
@@ -400,7 +405,9 @@ class MQuadrangle9 : public MQuadrangle {
   virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return QUA_9; }
+  virtual int getTypeForUNV(){ return 0; } // not available
   virtual char *getStringForPOS(){ return "SQ2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
 };
 
 class MTetrahedron : public MElement {
@@ -488,11 +495,17 @@ class MTetrahedron10 : public MTetrahedron {
     static const int map[10] = {0, 4, 1, 5, 2, 6, 8, 9, 7, 3};
     return getVertex(map[num]); 
   }
+  virtual MVertex *getVertexBDF(int num)
+  {
+    static const int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 6; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return TET_10; }
   virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
   virtual char *getStringForPOS(){ return "SS2"; }
+  virtual char *getStringForBDF(){ return "CTETRA"; }
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -591,14 +604,22 @@ class MHexahedron20 : public MHexahedron {
   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};
+    static const 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 MVertex *getVertexBDF(int num)
+  {
+    static const int map[20] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13, 
+				9, 10, 12, 14, 15, 16, 18, 19, 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 char *getStringForPOS(){ return 0; } // not available
+  virtual char *getStringForBDF(){ return "CHEXA"; }
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -646,7 +667,9 @@ class MHexahedron27 : public MHexahedron {
   virtual int getNumVolumeVertices(){ return 1; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return HEX_27; }
+  virtual int getTypeForUNV(){ return 0; } // not available
   virtual char *getStringForPOS(){ return "SH2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -753,13 +776,20 @@ class MPrism15 : public MPrism {
   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};
+    static const int map[15] = {0, 6, 1, 9, 2, 7, 8, 10, 11, 3, 12, 4, 14, 5, 13};
+    return getVertex(map[num]); 
+  }
+  virtual MVertex *getVertexBDF(int num)
+  {
+    static const int map[15] = {0, 1, 2, 3, 4, 5, 6, 9, 7, 8, 10, 11, 12, 14, 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 char *getStringForPOS(){ return 0; } // not available
+  virtual char *getStringForBDF(){ return "CPENTA"; }
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -800,7 +830,9 @@ class MPrism18 : public MPrism {
   virtual int getNumFaceVertices(){ return 3; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return PRI_18; }
+  virtual int getTypeForUNV(){ return 0; } // not available
   virtual char *getStringForPOS(){ return "SI2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -852,6 +884,7 @@ class MPyramid : public MElement {
 		   _v[quadfaces_pyramid[num - 4][3]]);
   }
   virtual int getTypeForMSH(){ return PYR_5; }
+  virtual int getTypeForUNV(){ return 0; } // not available
   virtual char *getStringForPOS(){ return "SY"; }
   virtual char *getStringForBDF(){ return "CPYRAM"; }
   virtual int getVolumeSign()
@@ -901,6 +934,9 @@ class MPyramid13 : public MPyramid {
   virtual int getNumEdgeVertices(){ return 8; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return PYR_13; }
+  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual char *getStringForPOS(){ return 0; } // not available
+  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
@@ -939,7 +975,9 @@ class MPyramid14 : public MPyramid {
   virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
   virtual int getTypeForMSH(){ return PYR_14; }
+  virtual int getTypeForUNV(){ return 0; } // not available
   virtual char *getStringForPOS(){ return "SY2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
   virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){