diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index a5f26b76741cf0cca8afe882c42dfe38a891df1e..737de44f56e0a6b0092328778f81c70f86bfdc09 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.78 2006-08-19 18:48:06 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.79 2006-09-03 07:44:09 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -80,7 +80,7 @@ void Print_Usage(char *name){
   Msg(DIRECT, "  -1, -2, -3            Perform 1D, 2D or 3D mesh generation, then exit");
   Msg(DIRECT, "  -saveall              Save all elements (discard physical group definitions)");
   Msg(DIRECT, "  -o file               Specify mesh output file name");
-  Msg(DIRECT, "  -format string        Set output mesh format (msh, unv, mesh, stl, vrml)");
+  Msg(DIRECT, "  -format string        Set output mesh format (msh, unv, bdf, mesh, stl, vrml)");
   Msg(DIRECT, "  -algo string          Select mesh algorithm (iso, tri, aniso, netgen, tetgen)");
   Msg(DIRECT, "  -smooth int           Set number of mesh smoothing steps");
   Msg(DIRECT, "  -optimize             Optimize quality of tetrahedral elements");
@@ -397,6 +397,8 @@ void Get_Options(int argc, char *argv[])
             CTX.mesh.format = FORMAT_MSH;
           else if(!strcmp(argv[i], "unv"))
             CTX.mesh.format = FORMAT_UNV;
+          else if(!strcmp(argv[i], "bdf"))
+            CTX.mesh.format = FORMAT_BDF;
           else if(!strcmp(argv[i], "mesh"))
             CTX.mesh.format = FORMAT_MESH;
 	  else if(!strcmp(argv[i], "stl"))
diff --git a/Common/Context.h b/Common/Context.h
index 5d0c21d7297d31e73451aea83744b1c2e06260af..cc0edf3569305e85cfa2c7ab37f6a1b08d86603e 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -183,7 +183,7 @@ public :
       return val;
     }
     int oldxtrude, oldxtrude_recombine;
-    int allow_degenerated_extrude, save_all, stl_binary, msh_binary;
+    int allow_degenerated_extrude, save_all, stl_binary, msh_binary, bdf_field_format;
     char *triangle_options;
     int smooth_normals, reverse_all_normals;
     double angle_smooth_normals;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index fb23d0236701a1f3320051c1a78128d6eb81743f..b7fda784a228ecd5bc2215997805ee750c890d70 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -848,6 +848,8 @@ StringXNumber MeshOptions_Number[] = {
 
   { F|O, "BetaSmoothMetric" ,opt_mesh_beta_smooth_metric, 0.9 ,
     "Maximum ratio of two consecutive edge lengths" },
+  { F|O, "BdfFieldFormat" , opt_mesh_bdf_field_format , 1. , 
+    "Field format for Nastran BDF files (0=free, 1=small, 2=large)" },
 
   { F|O, "CharacteristicLengthFactor" , opt_mesh_lc_factor , 1.0 ,
     "Factor applied to all characteristic lengths (and background meshes)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 2c1fcd26a482ecebdcee3791e9bf542a93b9a889..d6272aa239f9029a8083c44be123c6e7c47cbdcf 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.308 2006-08-27 23:10:35 geuzaine Exp $
+// $Id: Options.cpp,v 1.309 2006-09-03 07:44:09 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -4669,6 +4669,17 @@ double opt_mesh_stl_binary(OPT_ARGS_NUM)
   return CTX.mesh.stl_binary;
 }
 
+double opt_mesh_bdf_field_format(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET){
+    CTX.mesh.bdf_field_format = (int)val;
+    if(CTX.mesh.bdf_field_format < 0 ||
+       CTX.mesh.bdf_field_format > 2)
+      CTX.mesh.bdf_field_format = 1;
+  }
+  return CTX.mesh.bdf_field_format;
+}
+
 double opt_mesh_nb_smoothing(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 49e099cc0179793a27ab2c1dca80e684758b16a6..5ae072b19e24b293a75697e934fe391676bcaf61 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -448,6 +448,7 @@ double opt_mesh_format(OPT_ARGS_NUM);
 double opt_mesh_msh_file_version(OPT_ARGS_NUM);
 double opt_mesh_msh_binary(OPT_ARGS_NUM);
 double opt_mesh_stl_binary(OPT_ARGS_NUM);
+double opt_mesh_bdf_field_format(OPT_ARGS_NUM);
 double opt_mesh_nb_smoothing(OPT_ARGS_NUM);
 double opt_mesh_stl_distance_tol(OPT_ARGS_NUM);
 double opt_mesh_nb_elem_per_rc(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index e0589192a932c3b3abfa6d4a26da3feadb759d3e..55463e2c1bbd9dfa8e187df909854aece3455bdb 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.459 2006-09-02 22:24:23 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.460 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -679,8 +679,7 @@ int _save_mesh(char *name)
 
 int _save_bdf(char *name)
 {
-  CreateOutputFile(name, FORMAT_BDF);
-  return 1;
+  return bdf_dialog(name);
 }
 
 int _save_vrml(char *name)
@@ -751,6 +750,7 @@ int _save_auto(char *name)
   case FORMAT_OPT     : return _save_options(name);
   case FORMAT_MSH     : return _save_msh(name);
   case FORMAT_UNV     : return _save_unv(name);
+  case FORMAT_BDF     : return _save_bdf(name);
   case FORMAT_STL     : return _save_stl(name);
   case FORMAT_PS      : return _save_ps(name);
   case FORMAT_EPS     : return _save_eps(name);
diff --git a/Fltk/GUI_Extras.cpp b/Fltk/GUI_Extras.cpp
index b55e30c392a3b791c2b9023f28bdd1bd73fd9418..876964156d6c018e52a4db7457a48164d195ba4d 100644
--- a/Fltk/GUI_Extras.cpp
+++ b/Fltk/GUI_Extras.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI_Extras.cpp,v 1.23 2006-09-02 22:24:23 geuzaine Exp $
+// $Id: GUI_Extras.cpp,v 1.24 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -731,6 +731,67 @@ int unv_dialog(char *name)
   return 0;
 }
 
+// save bdf dialog
+
+int bdf_dialog(char *name)
+{
+  struct _bdf_dialog{
+    Fl_Window *window;
+    Fl_Choice *c;
+    Fl_Button *ok, *cancel;
+  };
+  static _bdf_dialog *dialog = NULL;
+
+  static Fl_Menu_Item formatmenu[] = {
+    {"Free field", 0, 0, 0},
+    {"Small field", 0, 0, 0},
+    {"Long field", 0, 0, 0},
+    {0}
+  };
+
+  const int BH = 2 * CTX.fontsize + 1;
+  const int BB = 7 * CTX.fontsize;
+  const int WB = 7;
+
+  if(!dialog){
+    dialog = new _bdf_dialog;
+    int h = 3 * WB + 2 * BH, w = 2 * BB + 3 * WB, y = WB;
+    // not a "Dialog_Window" since it is modal 
+    dialog->window = new Fl_Double_Window(w, h, "BDF Options");
+    dialog->window->box(GMSH_WINDOW_BOX);
+    dialog->c = new Fl_Choice(WB, y, BB + BB / 2, BH, "Format"); y += BH;
+    dialog->c->menu(formatmenu);
+    dialog->c->align(FL_ALIGN_RIGHT);
+    dialog->ok = new Fl_Return_Button(WB, y + WB, BB, BH, "OK");
+    dialog->cancel = new Fl_Button(2 * WB + BB, y + WB, BB, BH, "Cancel");
+    dialog->window->set_modal();
+    dialog->window->end();
+    dialog->window->hotspot(dialog->window);
+  }
+  
+  dialog->c->value(CTX.mesh.bdf_field_format);
+  dialog->window->show();
+
+  while(dialog->window->shown()){
+    Fl::wait();
+    for (;;) {
+      Fl_Widget* o = Fl::readqueue();
+      if (!o) break;
+      if (o == dialog->ok) {
+	opt_mesh_bdf_field_format(0, GMSH_SET | GMSH_GUI, dialog->c->value());
+	CreateOutputFile(name, FORMAT_BDF);
+	dialog->window->hide();
+	return 1;
+      }
+      if (o == dialog->window || o == dialog->cancel){
+	dialog->window->hide();
+	return 0;
+      }
+    }
+  }
+  return 0;
+}
+
 // save stl dialog
 
 int stl_dialog(char *name)
diff --git a/Fltk/GUI_Extras.h b/Fltk/GUI_Extras.h
index 569e47da2a97149d462305b5e3ca1e290a4489a6..13bdab2c2f77c55f69a4a19f05be8df2ebd5b8cd 100644
--- a/Fltk/GUI_Extras.h
+++ b/Fltk/GUI_Extras.h
@@ -36,6 +36,7 @@ int gl2ps_dialog(char *filename, char *title, int format);
 int options_dialog(char *filename);
 int msh_dialog(char *filename);
 int unv_dialog(char *filename);
+int bdf_dialog(char *filename);
 int stl_dialog(char *filename);
 
 #endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index fb7c00dcba05318db2a0e96ad5e0291e09f209eb..9b577e92cf38a3532f112bf85c8a0ab1329833ba 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -140,9 +140,10 @@ class GModel
   int readMESH(const std::string &name);
   int writeMESH(const std::string &name, double scalingFactor=1.0);
 
-  // IO for Nastran Bulk Data File format (free field with comma separator)
+  // IO for Nastran Bulk Data File format
   int readBDF(const std::string &name);
-  int writeBDF(const std::string &name, double scalingFactor=1.0);
+  int writeBDF(const std::string &name, int format=0, bool saveAll=false, 
+	       double scalingFactor=1.0);
 
   // FIXME: this will be removed (and rewritten)
   smooth_normals *normals;
diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index 7e0ee8cd6a783f38bc479396a9ca752abc2da2fe..97dd342d2b42ace957e0e430fa180955c685154c 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.37 2006-09-02 22:24:24 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.38 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -173,6 +173,24 @@ static void storePhysicalTagsInEntities(GModel *m, int dim,
   }
 }
 
+static bool checkVertexIndex(int index, std::map<int, MVertex*> &map)
+{
+  if(!map.count(index)){
+    Msg(GERROR, "Wrong vertex index %d", index);
+    return false;
+  }
+  return true;
+}
+
+static bool checkVertexIndex(int index, std::vector<MVertex*> &vec)
+{
+  if(index < 0 || index > (int)(vec.size() - 1)){
+    Msg(GERROR, "Wrong vertex index %d", index);
+    return false;
+  }
+  return true;
+}
+
 static int getNumVerticesForElementTypeMSH(int type)
 {
   switch (type) {
@@ -432,12 +450,18 @@ int GModel::readMSH(const std::string &name)
 	  }
 	  int indices[30];
 	  for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]);
-	  if(vertexVector.size())
+	  if(vertexVector.size()){
+	    for(int j = 0; j < numVertices; j++)
+	      if(!checkVertexIndex(indices[j], vertexVector)) return 0;
 	    createElementMSH(this, num, type, physical, elementary, partition, 
 			     indices, vertexVector, points, elements, physicals);
-	  else
+	  }
+	  else{
+	    for(int j = 0; j < numVertices; j++)
+	      if(!checkVertexIndex(indices[j], vertexMap)) return 0;
 	    createElementMSH(this, num, type, physical, elementary, partition, 
 			     indices, vertexMap, points, elements, physicals);
+	  }
 	  if(progress && (i % progress == progress - 1))
 	    Msg(PROGRESS, "Read %d elements", i + 1);
 	}
@@ -462,12 +486,18 @@ int GModel::readMSH(const std::string &name)
 	    int elementary = (numTags > 1) ? data[4 - numTags + 1] : 0;
 	    int partition = (numTags > 2) ? data[4 - numTags + 2] : 0;
 	    int *indices = &data[numTags + 1];
-	    if(vertexVector.size())
+	    if(vertexVector.size()){
+	      for(int j = 0; j < numVertices; j++)
+		if(!checkVertexIndex(indices[j], vertexVector)) return 0;
 	      createElementMSH(this, num, type, physical, elementary, partition,
 			       indices, vertexVector, points, elements, physicals);
-	    else
+	    }
+	    else{
+	      for(int j = 0; j < numVertices; j++)
+		if(!checkVertexIndex(indices[j], vertexMap)) return 0;
 	      createElementMSH(this, num, type, physical, elementary, partition, 
 			       indices, vertexMap, points, elements, physicals);
+	    }
 	    if(progress && ((numElementsPartial + i) % progress == progress - 1))
 	      Msg(PROGRESS, "Read %d elements", i + 1);
 	  }
@@ -946,12 +976,9 @@ static int readElementsVRML(FILE *fp, std::map<int, MVertex*> &v, int region,
       idx.push_back(i);
     }
     else{
-      for(unsigned int j = 0; j < idx.size(); j++){
-	if(!v.count(idx[j])){
-	  Msg(GERROR, "Bad vertex index in VRML file");
-	  return 0;
-	}
-      }
+      for(unsigned int j = 0; j < idx.size(); j++)
+	if(!checkVertexIndex(idx[j], v)) return 0;
+      
       if(idx.size() < 2){
 	Msg(INFO, "Skipping %d-vertex element", (int)idx.size());
       }
@@ -1155,7 +1182,10 @@ int GModel::readUNV(const std::string &name)
 	    if(sscanf(buffer, "%d %d %d", &dum, &dum, &dum) != 3) break;
 	  }
 	  int n[30];
-	  for(int i = 0; i < numNodes; i++) if(!fscanf(fp, "%d", &n[i])) break;
+	  for(int i = 0; i < numNodes; i++){
+	    if(!fscanf(fp, "%d", &n[i])) return 0;
+	    if(!checkVertexIndex(n[i], vertices)) return 0;
+	  }
 	  int type_msh = 0;
 	  switch(type){
 	  case 11: case 21: case 22: case 31:                   type_msh = LGN1; break;
@@ -1274,9 +1304,11 @@ int GModel::readMESH(const std::string &name)
     return 0;
   }
 
-  std::map<int, MVertex*> vertices;
+  std::vector<MVertex*> vertices;
   std::map<int, std::vector<MElement*> > elements[3];
 
+  int nbv, nbe, n[30], cl;
+
   while(!feof(fp)) {
     if(!fgets(buffer, sizeof(buffer), fp)) break;
     if(buffer[0] != '#'){ // skip comments
@@ -1286,54 +1318,55 @@ int GModel::readMESH(const std::string &name)
       }
       else if(!strcmp(str, "Vertices")){
 	if(!fgets(buffer, sizeof(buffer), fp)) break;
-	int nbv;
 	sscanf(buffer, "%d", &nbv);
 	Msg(INFO, "%d vertices", nbv);
+	vertices.resize(nbv);
 	for(int i = 0; i < nbv; i++) {
 	  if(!fgets(buffer, sizeof(buffer), fp)) break;
-	  int cl;
 	  double x, y, z;
 	  sscanf(buffer, "%lf %lf %lf %d", &x, &y, &z, &cl);
-	  vertices[i + 1] = new MVertex(x, y, z);
+	  vertices[i] = new MVertex(x, y, z);
 	}
       }
       else if(!strcmp(str, "Triangles")){
 	if(!fgets(buffer, sizeof(buffer), fp)) break;
-	int nbt;
-	sscanf(buffer, "%d", &nbt);
-	Msg(INFO, "%d triangles", nbt);
-	for(int i = 0; i < nbt; i++) {
+	sscanf(buffer, "%d", &nbe);
+	Msg(INFO, "%d triangles", nbe);
+	for(int i = 0; i < nbe; i++) {
 	  if(!fgets(buffer, sizeof(buffer), fp)) break;
-	  int n1, n2, n3, cl;
-	  sscanf(buffer, "%d %d %d %d", &n1, &n2, &n3, &cl);
+	  sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &cl);
+	  for(int j = 0; j < 3; j++)
+	    if(!checkVertexIndex(n[j] - 1, vertices)) return 0;
 	  elements[0][cl].push_back
-	    (new MTriangle(vertices[n1], vertices[n2], vertices[n3]));
+	    (new MTriangle(vertices[n[0] - 1], vertices[n[1] - 1], vertices[n[2] - 1]));
 	}
       }
       else if(!strcmp(str, "Quadrilaterals")) {
 	if(!fgets(buffer, sizeof(buffer), fp)) break;
-	int nbq;
-	sscanf(buffer, "%d", &nbq);
-	Msg(INFO, "%d quadrangles", nbq);
-	for(int i = 0; i < nbq; i++) {
+	sscanf(buffer, "%d", &nbe);
+	Msg(INFO, "%d quadrangles", nbe);
+	for(int i = 0; i < nbe; i++) {
 	  if(!fgets(buffer, sizeof(buffer), fp)) break;
-	  int n1, n2, n3, n4, cl;
-	  sscanf(buffer, "%d %d %d %d %d", &n1, &n2, &n3, &n4, &cl);
+	  sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
+	  for(int j = 0; j < 4; j++) 
+	    if(!checkVertexIndex(n[j] - 1, vertices)) return 0;
 	  elements[1][cl].push_back
-	    (new MQuadrangle(vertices[n1], vertices[n2], vertices[n3], vertices[n4]));
+	    (new MQuadrangle(vertices[n[0] - 1], vertices[n[1] - 1], 
+			     vertices[n[2] - 1], vertices[n[3] - 1]));
 	}
       }
       else if(!strcmp(str, "Tetrahedra")) {
 	if(!fgets(buffer, sizeof(buffer), fp)) break;
-	int nbt;
-	sscanf(buffer, "%d", &nbt);
-	Msg(INFO, "%d tetrahedra", nbt);
-	for(int i = 0; i < nbt; i++) {
+	sscanf(buffer, "%d", &nbe);
+	Msg(INFO, "%d tetrahedra", nbe);
+	for(int i = 0; i < nbe; i++) {
 	  if(!fgets(buffer, sizeof(buffer), fp)) break;
-	  int n1, n2, n3, n4, cl;
-	  sscanf(buffer, "%d %d %d %d %d", &n1, &n2, &n3, &n4, &cl);
+	  sscanf(buffer, "%d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &cl);
+	  for(int j = 0; j < 4; j++) 
+	    if(!checkVertexIndex(n[j] - 1, vertices)) return 0;
 	  elements[2][cl].push_back
-	    (new MTetrahedron(vertices[n1], vertices[n2], vertices[n3], vertices[n4]));
+	    (new MTetrahedron(vertices[n[0] - 1], vertices[n[1] - 1], 
+			      vertices[n[2] - 1], vertices[n[3] - 1]));
 	}
       }
     }
@@ -1412,38 +1445,108 @@ int GModel::writeMESH(const std::string &name, double scalingFactor)
   return 1;
 }
 
-static int readElementBDF(FILE *fp, char *buffer, unsigned int numNodes, 
-			  int *num, int *region, int *nodes)
+static int getFormatBDF(char *buffer, int keySize)
+{
+  if(buffer[keySize] == '*') return 2; // long fields
+  for(unsigned int i = 0; i < strlen(buffer); i++)
+    if(buffer[i] == ',') return 0; // free fields
+  return 1; // small fields;
+}
+
+static int readVertexBDF(FILE *fp, char *buffer, int keySize, 
+			 int *num, double *x, double *y, double *z)
+{
+  int dum;
+  char tmp[32], buffer2[256];
+
+  switch(getFormatBDF(buffer, keySize)){
+  case 0: // free field
+    if(sscanf(&buffer[keySize + 1], "%d,%d,%lf,%lf,%lf", num, &dum, x, y, z) != 5) 
+      return 0;
+    break;
+  case 1: // small field
+    tmp[8] = '\0';
+    strncpy(tmp, &buffer[8], 8); if(!sscanf(tmp, "%d", num)) return 0;
+    strncpy(tmp, &buffer[24], 8); if(!sscanf(tmp, "%lf", x)) return 0;
+    strncpy(tmp, &buffer[32], 8); if(!sscanf(tmp, "%lf", y)) return 0;
+    strncpy(tmp, &buffer[40], 8); if(!sscanf(tmp, "%lf", z)) return 0;
+    break;
+  case 2: // long field
+    tmp[16] = '\0';
+    strncpy(tmp, &buffer[8], 16); if(!sscanf(tmp, "%d", num)) return 0;
+    strncpy(tmp, &buffer[40], 16); if(!sscanf(tmp, "%lf", x)) return 0;
+    strncpy(tmp, &buffer[56], 16); if(!sscanf(tmp, "%lf", y)) return 0;
+    for(unsigned int i = 0; i < sizeof(buffer2); i++) buffer2[i] = '\0';
+    if(!fgets(buffer2, sizeof(buffer2), fp)) return 0;
+    // should check that continuation tags match
+    strncpy(tmp, &buffer2[8], 16); if(!sscanf(tmp, "%lf", z)) return 0;
+    break;
+  }
+  return 1;
+}
+
+static int readElementBDF(FILE *fp, char *buffer, int keySize, unsigned int numNodes, 
+			  int *num, int *region, std::vector<MVertex*> &vertices,
+			  std::map<int, MVertex*> &vertexMap)
 {
-  int newline = 0; 
+  // This routine is not general: when there is a continuation tag we
+  // just parse the next line. To follow the spec,
+  // 1. we should parse the whole remainder of the file to find
+  //    the corresponding continuation tag
+  // 2. we should parse more than one continuation line for elements
+  //    with more than 14 nodes (since we don't import such elements at
+  //    the moment parsing just one continuation line is OK)
+
+  char buffer2[256];
   std::vector<char*> vals;
+  int format = getFormatBDF(buffer, keySize);
+  int cmax = (format == 2) ? 16 : 8; // max char per (center) field
+  int nmax = (format == 2) ? 4 : 8; // max num of (center) fields per line
 
-  for(unsigned int i = 0; i < strlen(buffer); i++){
-    if(buffer[i] == ',') vals.push_back(&buffer[i + 1]);
-    else if(buffer[i] == '+'){ // the data continues on the next line
-      vals.pop_back();
-      newline = 1;
-      break;
+  for(unsigned int i = 0; i < sizeof(buffer2); i++) buffer2[i] = '\0';
+
+  if(format == 0){ // free fields
+    for(unsigned int i = 0; i < strlen(buffer); i++){
+      if(buffer[i] == ',') vals.push_back(&buffer[i + 1]);
+      if(vals.size() == 9) vals.pop_back(); // continuation field
+    }
+    if(vals.size() < 2 + numNodes){
+      if(!fgets(buffer2, sizeof(buffer2), fp)) return 0;
+      for(unsigned int i = 0; i < strlen(buffer2); i++){
+	if(buffer2[i] == ',') vals.push_back(&buffer2[i + 1]);
+      }
     }
   }
-  
-  if(newline){
-    char buffer2[256];
-    if(!fgets(buffer2, sizeof(buffer2), fp)) return 0;
-    for(unsigned int i = 0; i < strlen(buffer2); i++){
-      if(buffer2[i] == ',') vals.push_back(&buffer2[i + 1]);
+  else{ // small or long fields
+    for(int i = 0; i < nmax; i++){
+      if(vals.size() < 2 + numNodes) vals.push_back(&buffer[8 + cmax * i]);
     }
+    if(vals.size() < 2 + numNodes){
+      if(!fgets(buffer2, sizeof(buffer2), fp)) return 0;
+      for(int i = 0; i < nmax; i++){
+	if(vals.size() < 2 + numNodes) vals.push_back(&buffer2[8 + cmax * i]);
+      }
+    }    
   }
   
-  if(vals.size() < numNodes + 2){
+  if(vals.size() < 2 + numNodes){
     Msg(GERROR, "Missing nodes for element (%d < numNodes)", vals.size() - 2);
     return 0;
   }
-  
-  sscanf(vals[0], "%d", num);
-  sscanf(vals[1], "%d", region);
-  for(unsigned int i = 0; i < numNodes; i++)
-    if(!sscanf(vals[i + 2], "%d", &nodes[i])) return 0;
+
+  int n[30];
+  char tmp[32];
+  tmp[cmax] = '\0';
+  strncpy(tmp, vals[0], cmax); if(!sscanf(tmp, "%d", num)) *num = 0;
+  strncpy(tmp, vals[1], cmax); if(!sscanf(tmp, "%d", region)) *region = 0;
+  for(unsigned int i = 0; i < numNodes; i++){
+    strncpy(tmp, vals[i + 2], cmax); if(!sscanf(tmp, "%d", &n[i])) return 0;
+  }
+
+  for(unsigned int i = 0; i < numNodes; i++){
+    if(!checkVertexIndex(n[i], vertexMap)) return 0;
+    vertices.push_back(vertexMap[n[i]]);
+  }
   return 1;
 }
 
@@ -1456,76 +1559,63 @@ int GModel::readBDF(const std::string &name)
   }
 
   char buffer[256];
-  int num, dummy, region, n[30];
-  double x, y, z;
-  bool comma = false;
-
   std::map<int, MVertex*> vertices;
   std::map<int, std::vector<MElement*> > elements[6];
 
+  // nodes can be defined after elements, so parse the file twice
+
   while(!feof(fp)) {
+    for(unsigned int i = 0; i < sizeof(buffer); i++) buffer[i] = '\0';
     if(!fgets(buffer, sizeof(buffer), fp)) break;
     if(buffer[0] != '$'){ // skip comments
-      if(!comma){ // check that we have a free format file with comma separator
-	for(unsigned int i = 0; i < strlen(buffer); i++){
-	  if(buffer[i] == ','){ 
-	    comma = true;
-	    break; 
-	  }
-	}
-	if(!comma){
-	  Msg(GERROR, "BDF reader only accepts free field format, comma separated");
-	  break;
-	}
-      }
       if(!strncmp(buffer, "GRID", 4)){
-	sscanf(&buffer[5], "%d , %d , %lf, %lf , %lf", &num, &dummy, &x, &y, &z);
+	int num;
+	double x, y, z;
+	if(!readVertexBDF(fp, buffer, 4, &num, &x, &y, &z)) break;
 	vertices[num] = new MVertex(x, y, z);
       }
-      else if(!strncmp(buffer, "CBAR", 4)){
-	if(readElementBDF(fp, &buffer[4], 2, &num, &region, n))
-	  elements[0][region].push_back
-	    (new MLine(vertices[n[0]], vertices[n[1]], num));
+    }
+  }
+  Msg(INFO, "%d vertices", vertices.size());
+
+  rewind(fp);
+  while(!feof(fp)) {
+    for(unsigned int i = 0; i < sizeof(buffer); i++) buffer[i] = '\0';
+    if(!fgets(buffer, sizeof(buffer), fp)) break;
+    if(buffer[0] != '$'){ // skip comments
+      int num, region;
+      std::vector<MVertex*> n;
+      if(!strncmp(buffer, "CBAR", 4)){
+	if(readElementBDF(fp, buffer, 4, 2, &num, &region, n, vertices))
+	  elements[0][region].push_back(new MLine(n, num));
       }
       else if(!strncmp(buffer, "CTRIA3", 6)){
-	if(readElementBDF(fp, &buffer[6], 3, &num, &region, n))
-	  elements[1][region].push_back
-	    (new MTriangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], num));
+	if(readElementBDF(fp, buffer, 6, 3, &num, &region, n, vertices))
+	  elements[1][region].push_back(new MTriangle(n, num));
       }
       else if(!strncmp(buffer, "CTRIA6", 6)){
-	if(readElementBDF(fp, &buffer[6], 6, &num, &region, n))
-	  elements[1][region].push_back
-	    (new MTriangle2(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
-			    vertices[n[3]], vertices[n[4]], vertices[n[5]], num));
+	if(readElementBDF(fp, buffer, 6, 6, &num, &region, n, vertices))
+	  elements[1][region].push_back(new MTriangle2(n, num));
       }
       else if(!strncmp(buffer, "CQUAD4", 6)){
-	if(readElementBDF(fp, &buffer[6], 4, &num, &region, n))
-	  elements[2][region].push_back
-	    (new MQuadrangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
-			     vertices[n[3]], num));
+	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertices))
+	  elements[2][region].push_back(new MQuadrangle(n, num));
       }
       else if(!strncmp(buffer, "CTETRA", 6)){
-	if(readElementBDF(fp, &buffer[6], 4, &num, &region, n))
-	  elements[3][region].push_back
-	    (new MTetrahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
-			      vertices[n[3]], num));
+	if(readElementBDF(fp, buffer, 6, 4, &num, &region, n, vertices))
+	  elements[3][region].push_back(new MTetrahedron(n, num));
       }
       else if(!strncmp(buffer, "CHEXA", 5)){
-	if(readElementBDF(fp, &buffer[5], 8, &num, &region, n))
-	  elements[4][region].push_back
-	    (new MHexahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
-			     vertices[n[3]], vertices[n[4]], vertices[n[5]], 
-			     vertices[n[6]], vertices[n[7]], num));
+	if(readElementBDF(fp, buffer, 5, 8, &num, &region, n, vertices))
+	  elements[4][region].push_back(new MHexahedron(n, num));
       }
       else if(!strncmp(buffer, "CPENTA", 6)){
-	if(readElementBDF(fp, &buffer[6], 6, &num, &region, n))
-	  elements[5][region].push_back
-	    (new MPrism(vertices[n[0]], vertices[n[1]], vertices[n[2]], 
-			vertices[n[3]], vertices[n[4]], vertices[n[5]], num));
+	if(readElementBDF(fp, buffer, 6, 6, &num, &region, n, vertices))
+	  elements[5][region].push_back(new MPrism(n, num));
       }
     }
   }
-
+  
   for(int i = 0; i < (int)(sizeof(elements)/sizeof(elements[0])); i++) 
     storeElementsInEntities(this, elements[i]);
   associateEntityWithVertices(this);
@@ -1535,7 +1625,8 @@ int GModel::readBDF(const std::string &name)
   return 1;
 }
 
-int GModel::writeBDF(const std::string &name, double scalingFactor)
+int GModel::writeBDF(const std::string &name, int format, bool saveAll, 
+		     double scalingFactor)
 {
   FILE *fp = fopen(name.c_str(), "w");
   if(!fp){
@@ -1549,34 +1640,34 @@ int GModel::writeBDF(const std::string &name, double scalingFactor)
 
   for(viter it = firstVertex(); it != lastVertex(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeBDF(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
   for(fiter it = firstFace(); it != lastFace(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeBDF(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeBDF(fp, format, scalingFactor);
   
   for(eiter it = firstEdge(); it != lastEdge(); ++it){
     for(unsigned int i = 0; i < (*it)->lines.size(); i++)
-      (*it)->lines[i]->writeBDF(fp, (*it)->tag());
+      (*it)->lines[i]->writeBDF(fp, format, (*it)->tag());
   }
   for(fiter it = firstFace(); it != lastFace(); ++it){
     for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-      (*it)->triangles[i]->writeBDF(fp, (*it)->tag());
+      (*it)->triangles[i]->writeBDF(fp, format, (*it)->tag());
     for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-      (*it)->quadrangles[i]->writeBDF(fp, (*it)->tag());
+      (*it)->quadrangles[i]->writeBDF(fp, format, (*it)->tag());
   }
   for(riter it = firstRegion(); it != lastRegion(); ++it){
     for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
-      (*it)->tetrahedra[i]->writeBDF(fp, (*it)->tag());
+      (*it)->tetrahedra[i]->writeBDF(fp, format, (*it)->tag());
     for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
-      (*it)->hexahedra[i]->writeBDF(fp, (*it)->tag());
+      (*it)->hexahedra[i]->writeBDF(fp, format, (*it)->tag());
     for(unsigned int i = 0; i < (*it)->prisms.size(); i++)
-      (*it)->prisms[i]->writeBDF(fp, (*it)->tag());
+      (*it)->prisms[i]->writeBDF(fp, format, (*it)->tag());
   }
   
   fprintf(fp, "ENDDATA\n");
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 1a6d564489a81b4a6ebc4f3438f06d2d799f71d3..b143f604be23245613b5d0042d7cb95242dadf80 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.14 2006-09-02 22:24:24 geuzaine Exp $
+// $Id: MElement.cpp,v 1.15 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -280,16 +280,32 @@ void MElement::writeMESH(FILE *fp, int elementary)
   fprintf(fp, " %d\n", elementary);
 }
 
-void MElement::writeBDF(FILE *fp, int elementary)
+void MElement::writeBDF(FILE *fp, int format, int elementary)
 {
   char *str = getStringForBDF();
   if(str){
-    fprintf(fp, "%s,%d,%d", str, _num, elementary);
     int n = getNumVertices();
-    for(int i = 0; i < n; i++)
-      fprintf(fp, ",%d", getVertex(i)->getNum());
-    if(n == 2) // CBAR
-      fprintf(fp, ",0.,0.,0.");
-    fprintf(fp, "\n");
+    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(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, "%-8s%-8s%-8s", "0.", "0.", "0.");
+      fprintf(fp, "\n");
+    }
   }
 }
diff --git a/Geo/MElement.h b/Geo/MElement.h
index b32c665acb93bf46b5ea341ec8e342954a908b57..654195e0e727ba65aae3d49b6a2f5be118eefd58 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -137,7 +137,7 @@ class MElement
   virtual void writeVRML(FILE *fp);
   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 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;
@@ -160,16 +160,16 @@ class MLine : public MElement {
   }
   ~MLine(){}
   virtual int getDim(){ return 1; }
-  inline int getNumVertices(){ return 2; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 2; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 1; }
   virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
   virtual int getNumFaces(){ return 0; }
   virtual MFace getFace(int num){ throw; }
-  int getTypeForMSH(){ return LGN1; }
-  int getTypeForUNV(){ return 21; } // linear beam
-  char *getStringForPOS(){ return "SL"; }
-  char *getStringForBDF(){ return "CBAR"; }
+  virtual int getTypeForMSH(){ return LGN1; }
+  virtual int getTypeForUNV(){ return 21; } // linear beam
+  virtual char *getStringForPOS(){ return "SL"; }
+  virtual char *getStringForBDF(){ return "CBAR"; }
 };
 
 class MLine2 : public MLine {
@@ -187,12 +187,12 @@ class MLine2 : public MLine {
     _vs[0] = v[2];
   }
   ~MLine2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 3; }
-  inline MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
-  inline int getNumEdgeVertices(){ return 1; }
-  int getNumEdgesRep(){ return 2; }
-  MEdge getEdgeRep(int num)
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 3; }
+  virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
+  virtual int getNumEdgeVertices(){ return 1; }
+  virtual int getNumEdgesRep(){ return 2; }
+  virtual MEdge getEdgeRep(int num)
   { 
     static int edges_lin2[2][2] = {
       {0, 2}, {2, 1},
@@ -201,10 +201,10 @@ class MLine2 : public MLine {
     int i1 = edges_lin2[num][1];
     return MEdge(i0 < 2 ? _v[i0] : _vs[i0 - 2], i1 < 2 ? _v[i1] : _vs[i1 - 2]);
   }
-  int getTypeForMSH(){ return LGN2; }
-  int getTypeForUNV(){ return 24; } // parabolic beam
-  char *getStringForPOS(){ return "SL2"; }
-  char *getStringForBDF(){ return 0; }
+  virtual int getTypeForMSH(){ return LGN2; }
+  virtual int getTypeForUNV(){ return 24; } // parabolic beam
+  virtual char *getStringForPOS(){ return "SL2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
 };
 
 class MTriangle : public MElement {
@@ -223,9 +223,9 @@ class MTriangle : public MElement {
   }
   ~MTriangle(){}
   virtual int getDim(){ return 2; }
-  inline int getNumVertices(){ return 3; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
-  inline void revert () 
+  virtual int getNumVertices(){ return 3; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual void revert () 
     {
       MVertex *vv = _v[0];
       _v[0] = _v[1];
@@ -241,10 +241,10 @@ class MTriangle : public MElement {
   { 
     return MFace(_v[0], _v[1], _v[2]); 
   }
-  int getTypeForMSH(){ return TRI1; }
-  int getTypeForUNV(){ return 91; } // thin shell linear triangle
-  char *getStringForPOS(){ return "ST"; }
-  char *getStringForBDF(){ return "CTRIA3"; }
+  virtual int getTypeForMSH(){ return TRI1; }
+  virtual int getTypeForUNV(){ return 91; } // thin shell linear triangle
+  virtual char *getStringForPOS(){ return "ST"; }
+  virtual char *getStringForBDF(){ return "CTRIA3"; }
 };
 
 class MTriangle2 : public MTriangle {
@@ -263,12 +263,12 @@ class MTriangle2 : public MTriangle {
     for(int i = 0; i < 3; i++) _vs[i] = v[3 + i];
   }
   ~MTriangle2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 6; }
-  inline MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
-  inline int getNumEdgeVertices(){ return 3; }
-  int getNumEdgesRep(){ return 6; }
-  MEdge getEdgeRep(int num)
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 6; }
+  virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
+  virtual int getNumEdgeVertices(){ return 3; }
+  virtual int getNumEdgesRep(){ return 6; }
+  virtual MEdge getEdgeRep(int num)
   { 
     static int edges_tri2[6][2] = {
       {0, 3}, {3, 1},
@@ -279,8 +279,8 @@ class MTriangle2 : public MTriangle {
     int i1 = edges_tri2[num][1];
     return MEdge(i0 < 3 ? _v[i0] : _vs[i0 - 3], i1 < 3 ? _v[i1] : _vs[i1 - 3]);
   }
-  int getNumFacesRep(){ return 4; }
-  MFace getFaceRep(int num)
+  virtual int getNumFacesRep(){ return 4; }
+  virtual MFace getFaceRep(int num)
   { 
     static int trifaces_tri2[4][3] = {
       {0, 3, 5},
@@ -295,10 +295,10 @@ class MTriangle2 : public MTriangle {
 		 i1 < 3 ? _v[i1] : _vs[i1 - 3],
 		 i2 < 3 ? _v[i2] : _vs[i2 - 3]);
   }
-  int getTypeForMSH(){ return TRI2; }
-  int getTypeForUNV(){ return 92; } // thin shell parabolic triangle
-  char *getStringForPOS(){ return "ST2"; }
-  char *getStringForBDF(){ return "CTRIA6"; }
+  virtual int getTypeForMSH(){ return TRI2; }
+  virtual int getTypeForUNV(){ return 92; } // thin shell parabolic triangle
+  virtual char *getStringForPOS(){ return "ST2"; }
+  virtual char *getStringForBDF(){ return "CTRIA6"; }
 };
 
 class MQuadrangle : public MElement {
@@ -317,8 +317,8 @@ class MQuadrangle : public MElement {
   }
   ~MQuadrangle(){}
   virtual int getDim(){ return 2; }
-  inline int getNumVertices(){ return 4; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 4; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 4; }
   virtual MEdge getEdge(int num)
   {
@@ -326,10 +326,10 @@ class MQuadrangle : public MElement {
   }
   virtual int getNumFaces(){ return 1; }
   virtual MFace getFace(int num){ return MFace(_v[0], _v[1], _v[2], _v[3]); }
-  int getTypeForMSH(){ return QUA1; }
-  int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
-  char *getStringForPOS(){ return "SQ"; }
-  char *getStringForBDF(){ return "CQUAD4"; }
+  virtual int getTypeForMSH(){ return QUA1; }
+  virtual int getTypeForUNV(){ return 94; } // thin shell linear quadrilateral
+  virtual char *getStringForPOS(){ return "SQ"; }
+  virtual char *getStringForBDF(){ return "CQUAD4"; }
 };
 
 class MQuadrangle2 : public MQuadrangle {
@@ -340,7 +340,7 @@ class MQuadrangle2 : public MQuadrangle {
 	       MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8, int num=0, int part=0) 
     : MQuadrangle(v0, v1, v2, v3, num, part)
   {
-    _vs[0] = v4; _v[1] = v5; _v[2] = v6; _v[3] = v7; _v[4] = v8;
+    _vs[0] = v4; _vs[1] = v5; _vs[2] = v6; _vs[3] = v7; _vs[4] = v8;
   }
   MQuadrangle2(std::vector<MVertex*> &v, int num=0, int part=0) 
     : MQuadrangle(v, num, part)
@@ -348,16 +348,16 @@ class MQuadrangle2 : public MQuadrangle {
     for(int i = 0; i < 5; i++) _vs[i] = v[4 + i];
   }
   ~MQuadrangle2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 9; }
-  inline MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
-  inline int getNumEdgeVertices(){ return 4; }
-  inline int getNumFaceVertices(){ return 1; }
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 9; }
+  virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual int getNumEdgeVertices(){ return 4; }
+  virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  int getTypeForMSH(){ return QUA2; }
-  int getTypeForUNV(){ return 95; } // thin shell parabolic quadrilateral
-  char *getStringForPOS(){ return "SQ2"; }
-  char *getStringForBDF(){ return 0; }
+  virtual int getTypeForMSH(){ return QUA2; }
+  virtual int getTypeForUNV(){ return 95; } // thin shell parabolic quadrilateral
+  virtual char *getStringForPOS(){ return "SQ2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
 };
 
 class MTetrahedron : public MElement {
@@ -376,8 +376,8 @@ class MTetrahedron : public MElement {
   }
   ~MTetrahedron(){}
   virtual int getDim(){ return 3; }
-  inline int getNumVertices(){ return 4; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 4; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 6; }
   virtual MEdge getEdge(int num)
   {
@@ -390,10 +390,10 @@ class MTetrahedron : public MElement {
 		 _v[trifaces_tetra[num][1]],
 		 _v[trifaces_tetra[num][2]]);
   }
-  int getTypeForMSH(){ return TET1; }
-  int getTypeForUNV(){ return 111; } // solid linear tetrahedron
-  char *getStringForPOS(){ return "SS"; }
-  char *getStringForBDF(){ return "CTETRA"; }
+  virtual int getTypeForMSH(){ return TET1; }
+  virtual int getTypeForUNV(){ return 111; } // solid linear tetrahedron
+  virtual char *getStringForPOS(){ return "SS"; }
+  virtual char *getStringForBDF(){ return "CTETRA"; }
   virtual double getVolume()
   { 
     double mat[3][3];
@@ -409,7 +409,7 @@ class MTetrahedron : public MElement {
     return det3x3(mat) / 6.;
   }
   virtual int getVolumeSign(){ return sign(getVolume()); }
-  void setVolumePositive()
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -437,16 +437,16 @@ class MTetrahedron2 : public MTetrahedron {
     for(int i = 0; i < 6; i++) _vs[i] = v[4 + i];
   }
   ~MTetrahedron2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 10; }
-  inline MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
-  inline int getNumEdgeVertices(){ return 6; }
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 10; }
+  virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+  virtual int getNumEdgeVertices(){ return 6; }
   // TODO: edgeRep, faceRep
-  int getTypeForMSH(){ return TET2; }
-  int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
-  char *getStringForPOS(){ return "SS2"; }
-  char *getStringForBDF(){ return 0; }
-  void setVolumePositive()
+  virtual int getTypeForMSH(){ return TET2; }
+  virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
+  virtual char *getStringForPOS(){ return "SS2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -475,8 +475,8 @@ class MHexahedron : public MElement {
   }
   ~MHexahedron(){}
   virtual int getDim(){ return 3; }
-  inline int getNumVertices(){ return 8; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 8; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 12; }
   virtual MEdge getEdge(int num)
   {
@@ -490,10 +490,10 @@ class MHexahedron : public MElement {
 		 _v[quadfaces_hexa[num][2]],
 		 _v[quadfaces_hexa[num][3]]);
   }
-  int getTypeForMSH(){ return HEX1; }
-  int getTypeForUNV(){ return 115; } // solid linear brick
-  char *getStringForPOS(){ return "SH"; }
-  char *getStringForBDF(){ return "CHEXA"; }
+  virtual int getTypeForMSH(){ return HEX1; }
+  virtual int getTypeForUNV(){ return 115; } // solid linear brick
+  virtual char *getStringForPOS(){ return "SH"; }
+  virtual char *getStringForBDF(){ return "CHEXA"; }
   virtual int getVolumeSign()
   { 
     double mat[3][3];
@@ -508,7 +508,7 @@ class MHexahedron : public MElement {
     mat[2][2] = _v[4]->z() - _v[0]->z();
     return sign(det3x3(mat));
   }
-  void setVolumePositive()
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -541,18 +541,18 @@ class MHexahedron2 : public MHexahedron {
     for(int i = 0; i < 19; i++) _vs[i] = v[8 + i];
   }
   ~MHexahedron2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 27; }
-  inline MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
-  inline int getNumEdgeVertices(){ return 12; }
-  inline int getNumFaceVertices(){ return 6; }
-  inline int getNumVolumeVertices(){ return 1; }
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 27; }
+  virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+  virtual int getNumEdgeVertices(){ return 12; }
+  virtual int getNumFaceVertices(){ return 6; }
+  virtual int getNumVolumeVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  int getTypeForMSH(){ return HEX2; }
-  int getTypeForUNV(){ return 116; } // solid parabolic brick
-  char *getStringForPOS(){ return "SH2"; }
-  char *getStringForBDF(){ return 0; }
-  void setVolumePositive()
+  virtual int getTypeForMSH(){ return HEX2; }
+  virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
+  virtual char *getStringForPOS(){ return "SH2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -586,8 +586,8 @@ class MPrism : public MElement {
   }
   ~MPrism(){}
   virtual int getDim(){ return 3; }
-  inline int getNumVertices(){ return 6; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 6; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 9; }
   virtual MEdge getEdge(int num)
   {
@@ -606,10 +606,10 @@ class MPrism : public MElement {
 		   _v[quadfaces_prism[num - 2][2]],
 		   _v[quadfaces_prism[num - 2][3]]);
   }
-  int getTypeForMSH(){ return PRI1; }
-  int getTypeForUNV(){ return 112; } // solid linear wedge
-  char *getStringForPOS(){ return "SI"; }
-  char *getStringForBDF(){ return "CPENTA"; }
+  virtual int getTypeForMSH(){ return PRI1; }
+  virtual int getTypeForUNV(){ return 112; } // solid linear wedge
+  virtual char *getStringForPOS(){ return "SI"; }
+  virtual char *getStringForBDF(){ return "CPENTA"; }
   virtual int getVolumeSign()
   { 
     double mat[3][3];
@@ -624,7 +624,7 @@ class MPrism : public MElement {
     mat[2][2] = _v[3]->z() - _v[0]->z();
     return sign(det3x3(mat));
   }
-  void setVolumePositive()
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -654,17 +654,17 @@ class MPrism2 : public MPrism {
     for(int i = 0; i < 12; i++) _vs[i] = v[6 + i];
   }
   ~MPrism2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 18; }
-  inline MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
-  inline int getNumEdgeVertices(){ return 9; }
-  inline int getNumFaceVertices(){ return 3; }
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 18; }
+  virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
+  virtual int getNumEdgeVertices(){ return 9; }
+  virtual int getNumFaceVertices(){ return 3; }
   // TODO: edgeRep, faceRep
-  int getTypeForMSH(){ return PRI2; }
-  int getTypeForUNV(){ return 113; } // solid parabolic wedge
-  char *getStringForPOS(){ return "SI2"; }
-  char *getStringForBDF(){ return 0; }
-  void setVolumePositive()
+  virtual int getTypeForMSH(){ return PRI2; }
+  virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
+  virtual char *getStringForPOS(){ return "SI2"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -694,8 +694,8 @@ class MPyramid : public MElement {
   }
   ~MPyramid(){}
   virtual int getDim(){ return 3; }
-  inline int getNumVertices(){ return 5; }
-  inline MVertex *getVertex(int num){ return _v[num]; }
+  virtual int getNumVertices(){ return 5; }
+  virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual int getNumEdges(){ return 8; }
   virtual MEdge getEdge(int num)
   {
@@ -714,10 +714,10 @@ class MPyramid : public MElement {
 		   _v[quadfaces_pyramid[num - 4][2]],
 		   _v[quadfaces_pyramid[num - 4][3]]);
   }
-  int getTypeForMSH(){ return PYR1; }
-  int getTypeForUNV(){ return 0; }
-  char *getStringForPOS(){ return "SY"; }
-  char *getStringForBDF(){ return 0; }
+  virtual int getTypeForMSH(){ return PYR1; }
+  virtual int getTypeForUNV(){ return 0; } // not available
+  virtual char *getStringForPOS(){ return "SY"; }
+  virtual char *getStringForBDF(){ return 0; } // not available
   virtual int getVolumeSign()
   { 
     double mat[3][3];
@@ -732,7 +732,7 @@ class MPyramid : public MElement {
     mat[2][2] = _v[4]->z() - _v[0]->z();
     return sign(det3x3(mat));
   }
-  void setVolumePositive()
+  virtual void setVolumePositive()
   {
     if(getVolumeSign() < 0){
       MVertex *tmp;
@@ -760,17 +760,17 @@ class MPyramid2 : public MPyramid {
     for(int i = 0; i < 9; i++) _vs[i] = v[5 + i];
   }
   ~MPyramid2(){}
-  inline int getPolynomialOrder(){ return 2; }
-  inline int getNumVertices(){ return 14; }
-  inline MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
-  inline int getNumEdgeVertices(){ return 8; }
-  inline int getNumFaceVertices(){ return 1; }
+  virtual int getPolynomialOrder(){ return 2; }
+  virtual int getNumVertices(){ return 14; }
+  virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
+  virtual int getNumEdgeVertices(){ return 8; }
+  virtual int getNumFaceVertices(){ return 1; }
   // TODO: edgeRep, faceRep
-  int getTypeForMSH(){ return PYR2; }
-  int getTypeForUNV(){ return 0; }
-  char *getStringForPOS(){ return "SY2"; }
-  char *getStringForBDF(){ return 0; }
-  void setVolumePositive()
+  virtual int getTypeForMSH(){ return PYR2; }
+  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){
       MVertex *tmp;
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 2875fff59882da0576cc89b3bd50fdac9251e81b..369c395d065114cb1fb417653b11d96d7ee1490a 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.8 2006-09-02 22:24:24 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.9 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -83,12 +83,45 @@ void MVertex::writeUNV(FILE *fp, double scalingFactor)
 
 void MVertex::writeMESH(FILE *fp, double scalingFactor)
 {
-  fprintf(fp, " %20.14E      %20.14E      %20.14E      %d\n", 
+  fprintf(fp, " %20.14G      %20.14G      %20.14G      %d\n", 
 	  x() * scalingFactor, y() * scalingFactor, z() * scalingFactor, 0);
 }
 
-void MVertex::writeBDF(FILE *fp, double scalingFactor)
+static void double_to_char8(double val, char *str){
+  if(val >= 1.e6)
+    sprintf(str, "%.2E", val);
+  else if(val >= 1.e-3)
+    sprintf(str, "%f", val);
+  else if(val >= 0)
+    sprintf(str, "%.2E", val);
+  else if(val >= -1.e-3)
+    sprintf(str, "%.1E", val);
+  else if(val >= -1.e6)
+    sprintf(str, "%f", val);
+  else
+    sprintf(str, "%.1E", val);
+  str[8] = '\0';
+}
+
+void MVertex::writeBDF(FILE *fp, int format, double scalingFactor)
 {
-  fprintf(fp, "GRID,%d,%d,%f,%f,%f\n", _num, 0,
-	  x() * scalingFactor, y() * scalingFactor, z() * scalingFactor);
+  char xs[17], ys[17], zs[17];
+  double x1 = x() * scalingFactor;
+  double y1 = y() * scalingFactor;
+  double z1 = z() * scalingFactor;
+  if(format == 0){
+    // free field format (max 8 char per field, comma separated, 10 per line)
+    double_to_char8(x1, xs); double_to_char8(y1, ys); double_to_char8(z1, zs);
+    fprintf(fp, "GRID,%d,%d,%s,%s,%s\n", _num, 0, xs, ys, zs);
+  }
+  else if(format == 1){ 
+    // small field format (8 char par field, 10 per line)
+    double_to_char8(x1, xs); double_to_char8(y1, ys); double_to_char8(z1, zs);
+    fprintf(fp, "GRID    %-8d%-8d%-8s%-8s%-8s\n", _num, 0, xs, ys, zs);
+  }
+  else{ 
+    // large field format (8 char first/last field, 16 char middle, 6 per line)
+    fprintf(fp, "GRID*   %-16d%-16d%-16.9G%-16.9G*N%-6d\n", _num, 0, x1, y1, _num);
+    fprintf(fp, "*N%-6d%-16.9G\n", _num, z1);
+  }
 }
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 688c8e82bb484c4d64c27b2d14d198d3677b20ee..c6f32f211a3dfaa2ddcf612b98e17005564bdd28 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -83,7 +83,7 @@ class MVertex{
   void writeVRML(FILE *fp, double scalingFactor=1.0);
   void writeUNV(FILE *fp, double scalingFactor=1.0);
   void writeMESH(FILE *fp, double scalingFactor=1.0);
-  void writeBDF(FILE *fp, double scalingFactor=1.0);
+  void writeBDF(FILE *fp, int format=0, double scalingFactor=1.0);
 };
 
 class MEdgeVertex : public MVertex{
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index c6ae1946dc30ad56ec4fe8152b624e03f209b939..c1d3823cf33fca4cb46bfab0469987d028dea885 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.95 2006-09-01 10:10:05 remacle Exp $
+// $Id: Generator.cpp,v 1.96 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -399,8 +399,7 @@ void mai3d(int ask)
 
   CTX.threads_lock = 1;
 
-  // Clean up all the 2nd order nodes and transfer all SimplexBase
-  // into "real" Simplexes
+  // Change any high order elements back into first order ones
   Degre1();
 
   // 1D mesh
diff --git a/Mesh/Makefile b/Mesh/Makefile
index f8a30e427ab4a179f9a1948bd5790b7c024e3e65..76fe131e9839459e579602b360317ed1ee2cd2d5 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.127 2006-09-01 10:10:05 remacle Exp $
+# $Id: Makefile,v 1.128 2006-09-03 07:44:10 geuzaine Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -360,28 +360,16 @@ Generator.o: Generator.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../Common/ColorTable.h ../Common/VertexArray.h \
   ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
   ../Common/GmshMatrix.h Create.h ../Parser/OpenFile.h PartitionMesh.h \
-  ../Common/OS.h meshGEdge.h meshGFace.h ../Geo/GModel.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/MVertex.h ../Geo/GPoint.h ../Geo/GEdge.h \
-  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/GFace.h ../Geo/GRegion.h \
-  ../Geo/GEntity.h ../Geo/MElement.h ../Geo/SBoundingBox3d.h
-# 1 "/Users/geuzaine/.gmsh/Mesh//"
-DiscreteSurface.o: DiscreteSurface.cpp ../Common/Gmsh.h \
-  ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../Numeric/Numeric.h Mesh.h ../Common/GmshDefines.h \
-  Vertex.h Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h \
-  Metric.h Matrix.h ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Geo/ExtrudeParams.h ../Geo/Geo.h Create.h Interpolation.h \
-  ../Common/Context.h BDS.h ../Geo/GFace.h ../Geo/GPoint.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/MElement.h \
-  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
-  ../Common/GmshMatrix.h PartitionMesh.h ../Common/OS.h
+  ../Common/OS.h meshGEdge.h meshGFace.h meshGRegion.h ../Geo/GModel.h \
+  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/MVertex.h ../Geo/GPoint.h \
+  ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/GFace.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \
+  ../Geo/SBoundingBox3d.h
+# 1 "/Users/geuzaine/.gmsh/Mesh//"
+DiscreteSurface.o: DiscreteSurface.cpp Mesh.h ../Common/GmshDefines.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h Vertex.h Element.h \
+  Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h Metric.h Matrix.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
 SwapEdge.o: SwapEdge.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -429,23 +417,35 @@ meshGEdge.o: meshGEdge.cpp meshGEdge.h ../Geo/GEdge.h ../Geo/GEntity.h \
   Utils.h Vertex.h Mesh.h Element.h Simplex.h Face.h Edge.h \
   ../Geo/ExtrudeParams.h Metric.h Matrix.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
-meshGFace.o: meshGFace.cpp meshGFace.h ../Geo/SPoint2.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
-  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GPoint.h 2D_Mesh.h \
-  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h Mesh.h Vertex.h \
-  Element.h Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h Metric.h \
-  Matrix.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/MElement.h ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
-  ../Numeric/Numeric.h ../Common/Context.h ../Geo/GFace.h ../Geo/GPoint.h \
-  ../Geo/GEntity.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h Utils.h ../Common/Message.h BDS.h ../Common/Views.h \
+meshGFace.o: meshGFace.cpp meshGFace.h ../Geo/GVertex.h ../Geo/GEntity.h \
+  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+  ../Geo/SPoint3.h ../Common/GmshDefines.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h 2D_Mesh.h ../DataStr/List.h \
+  ../DataStr/Tree.h ../DataStr/avl.h Mesh.h Vertex.h Element.h Simplex.h \
+  Face.h Edge.h ../Geo/ExtrudeParams.h Metric.h Matrix.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/MVertex.h \
+  ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/Numeric.h \
+  ../Common/Context.h ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h \
+  ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  Utils.h ../Common/Message.h BDS.h ../Common/Views.h \
   ../Common/ColorTable.h ../Common/VertexArray.h \
   ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
   ../Common/GmshMatrix.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
+meshGRegion.o: meshGRegion.cpp meshGRegion.h ../Geo/GRegion.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
+  ../Geo/MElement.h ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
+  ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/Numeric.h \
+  ../Common/Context.h ../DataStr/List.h ../Geo/GFace.h ../Geo/GPoint.h \
+  ../Geo/GEntity.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h BDS.h ../Common/Views.h ../Common/ColorTable.h \
+  ../Common/VertexArray.h ../Common/SmoothNormals.h \
+  ../Common/AdaptiveViews.h ../Common/GmshMatrix.h ../Common/Message.h
+# 1 "/Users/geuzaine/.gmsh/Mesh//"
 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 ../DataStr/List.h ../DataStr/Tree.h \
@@ -478,22 +478,10 @@ SecondOrder.o: SecondOrder.cpp ../Geo/GModel.h ../Geo/GVertex.h \
   ../Geo/MEdge.h ../Geo/MElement.h ../Common/VertexArray.h \
   ../Common/Message.h ../Common/OS.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
-PartitionMesh.o: PartitionMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
-  ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
-  ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
-  ../Numeric/Numeric.h Mesh.h ../Common/GmshDefines.h Vertex.h Element.h \
+PartitionMesh.o: PartitionMesh.cpp Mesh.h ../Common/GmshDefines.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h Vertex.h Element.h \
   Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h Metric.h Matrix.h \
-  ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Geo/ExtrudeParams.h \
-  ../Geo/Geo.h Create.h Interpolation.h ../Common/Context.h BDS.h \
-  ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Geo/MElement.h ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h ../Common/Views.h ../Common/ColorTable.h \
-  ../Common/VertexArray.h ../Common/SmoothNormals.h \
-  ../Common/AdaptiveViews.h ../Common/GmshMatrix.h PartitionMesh.h \
-  ../Parser/OpenFile.h
+  PartitionMesh.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
 Smoothing.o: Smoothing.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Parser/CreateFile.cpp b/Parser/CreateFile.cpp
index 9bd386e43ed18a1a9f5f276cb72ca4a3693aeb38..39996d26dfe9d3c62aa2588be6ae0f32aa463a81 100644
--- a/Parser/CreateFile.cpp
+++ b/Parser/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.7 2006-09-02 22:24:24 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.8 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -61,6 +61,7 @@ int GuessFileFormatFromFileName(char *name)
   else if(!strcmp(ext, ".stl"))     return FORMAT_STL;
   else if(!strcmp(ext, ".mesh"))    return FORMAT_MESH;
   else if(!strcmp(ext, ".bdf"))     return FORMAT_BDF;
+  else if(!strcmp(ext, ".nas"))     return FORMAT_BDF;
   else if(!strcmp(ext, ".wrl"))     return FORMAT_VRML;
   else if(!strcmp(ext, ".vrml"))    return FORMAT_VRML;
   else if(!strcmp(ext, ".gif"))     return FORMAT_GIF;
@@ -166,7 +167,8 @@ void CreateOutputFile(char *filename, int format)
     break;
 
   case FORMAT_BDF:
-    GMODEL->writeBDF(name, CTX.mesh.scaling_factor);
+    GMODEL->writeBDF(name, CTX.mesh.bdf_field_format, CTX.mesh.save_all,
+		     CTX.mesh.scaling_factor);
     break;
 
   case FORMAT_POS:
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index a239547d42192de7ffd34ba5b5cd9eeae6b1ce97..a2b23531c12456c9aaddda58f119c73ad0f3842b 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.120 2006-09-02 22:24:24 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.121 2006-09-03 07:44:10 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -296,7 +296,8 @@ int MergeProblem(char *name, int warn_if_missing)
   else if(!strcmp(ext, ".mesh") || !strcmp(ext, ".MESH")){
     status = GMODEL->readMESH(name);
   }
-  else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF")){
+  else if(!strcmp(ext, ".bdf") || !strcmp(ext, ".BDF") ||
+	  !strcmp(ext, ".nas") || !strcmp(ext, ".NAS")){
     status = GMODEL->readBDF(name);
   }
 #if defined(HAVE_FLTK)
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 9ba8730efc4e6909c987a55c6dada7342bce6c65..67c21439f8b5095b6e5be577eff33c2a508a2d1c 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,10 +1,11 @@
-$Id: VERSIONS,v 1.364 2006-08-26 17:00:25 geuzaine Exp $
+$Id: VERSIONS,v 1.365 2006-09-03 07:44:11 geuzaine Exp $
 
 2.0: new geometry and mesh databases; complete rewrite of geometry and
-mesh drawing code; complete rewrite of the input/output code (new
-binary mesh format, improved support for STL, MESH, VRML and UNV
-formats); new 2D mesh algorithm; option changes in the interface are
-now applied instantenously; lots of improvements all over the place.
+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.
 
 1.66: added support for offscreen rendering using OSMesa; added
 support for SVG output;
diff --git a/doc/texinfo/command_line.texi b/doc/texinfo/command_line.texi
index 1d263678d42e1d36bfb67a43f9c7232833c90bd5..82886b8ae1a27eceadeadd7805bbad70539b6be2 100644
--- a/doc/texinfo/command_line.texi
+++ b/doc/texinfo/command_line.texi
@@ -17,7 +17,7 @@ Save all elements (discard physical group definitions)
 @item -o file
 Specify mesh output file name
 @item -format string
-Set output mesh format (msh, unv, mesh, stl, vrml)
+Set output mesh format (msh, unv, bdf, mesh, stl, vrml)
 @item -algo string
 Select mesh algorithm (iso, tri, aniso, netgen, tetgen)
 @item -smooth int
diff --git a/utils/converters/nastran/nastran2msh.c b/utils/converters/nastran/nastran2msh.c
deleted file mode 100644
index 051225ebc325f204ef32c26b2f7b884f67d7b057..0000000000000000000000000000000000000000
--- a/utils/converters/nastran/nastran2msh.c
+++ /dev/null
@@ -1,180 +0,0 @@
-/* Auteur: Jean Vis (---.magotteaux.com) */
-
-/*Include Linux:*/
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-
-#define VERSION "0.01: 01 dec 2005"
-
-
-struct noeuds {
- float x,y,z;
- long int no;
-} noeudsl[1000000];
-
-struct elem {
-  long int type,n,grp,n1,n2,n3,n4,n5,n6,n7,n8;
-} elemsl[1000000];
-
-long int nbN,nbE;
-double xf,yf,zf;
-double cfx,cfy,gl,frot;
-float xb=425.0,yb=0.0;
-float xm,ym;
-char basedata[256],outdata[300],version[17]=VERSION,enregistre[246];
-
-/*#############################################################*/
-/*#############################################################*/
-int datain()
-{
- FILE *ft;
- long int i,j,nonoe,ts,tts,noelem,ngrp,n1,n2,n3,n4,n5,n6,n7,n8;
- double n,x,y,z;
- char typeitm[16],comm[256],ligne[256],comm2[256],f3[30];
- char xtxt[20],ytxt[20],ztxt[20];
-
- printf("Fichier de donnees:");
- gets(basedata);
-
-
-/*Lecture fichier in */
-/*********************/
-
- printf("\nLecture des Noeuds: \n");
- ft=fopen(basedata,"r");
- if (ft==NULL) return (-1);
- nbN=0;
- nbE=0;
- while (!feof(ft))
- {
-   fgets(ligne,sizeof ligne,ft);
-   if (strncmp(ligne,"$",1)){
-      if (!strncmp(ligne,"GRID",4)){
-/* ici on a une ligne qui commence par "GRID" cad un noeud*/      
-      
-/*printf("ligne:%i ->%s",i,ligne);*/
-
-         ts=sscanf(ligne,"%s %ld %24s",comm,&nonoe,f3);
-/*printf("ligne:%i (ts:%i)->no=%i (%s) 3f=%s\n",i,ts,nonoe,comm,f3);*/
-         strncpy(xtxt,f3,7);
-         strncpy(ytxt,f3+8,8);
-         strncpy(ztxt,f3+16,8);
-/*printf ("x: %s y: %s z: %s\n",xtxt,ytxt,ztxt);*/
-         x=atof(xtxt);
-         y=atof(ytxt);
-         z=atof(ztxt);
-/*printf ("x: %f y: %f z: %f\n\n",x,y,z);*/
-         noeudsl[nbN].x=x;
-         noeudsl[nbN].y=y;
-         noeudsl[nbN].z=z;
-         noeudsl[nbN].no=nonoe;
-         nbN++;
-      }
-
-      if (!strncmp(ligne,"CTRIA3",6)){
-/*printf("ligne:%i ->%s",i,ligne);*/
-         ts=sscanf(ligne,"%s %d %d %d %d %d",comm,&noelem,&ngrp,&n1,&n2,&n3);
-         elemsl[nbE].type=2;
-         elemsl[nbE].n=noelem;
-         elemsl[nbE].grp=ngrp;
-         elemsl[nbE].n1=n1;
-         elemsl[nbE].n2=n2;
-         elemsl[nbE].n3=n3;
-         nbE++;
-      }
-
-      if (!strncmp(ligne,"CTETRA",6)){
-/*printf("ligne:%i ->%s",i,ligne);*/
-         ts=sscanf(ligne,"%s %d %d %d %d %d %d",comm,&noelem,&ngrp,&n1,&n2,&n3,&n4);
-         elemsl[nbE].type=4;
-         elemsl[nbE].n=noelem;
-         elemsl[nbE].grp=ngrp;
-         elemsl[nbE].n1=n1;
-         elemsl[nbE].n2=n2;
-         elemsl[nbE].n3=n3;
-         elemsl[nbE].n4=n4;
-         nbE++;
-      }
-      
-   }
-   i++;
- }
- 
- fclose(ft);
-
- return (0);
-}
-/*#############################################################*/
-
-/*#############################################################*/
-int dataout()
-{
-
-/*################################################
-  #                                              #
-  #  Goal: write a 1.0 Version .msh file of gmsh #
-  #                                              #
-  ################################################
-*/
-
- FILE *ft;
- long int i;
-
-/*Construction nom de fichier gmsh */
- strcpy(outdata,basedata);
- strcat(outdata,".msh");
-
- ft=fopen(outdata,"w");
- 
-/* Ecriture des noeuds 
-######################*/
-
- fprintf (ft,"$NOD\n");
- fprintf (ft,"%ld\n",nbN);
-
- i=0;
- while (i <nbN)
- {
-    fprintf (ft,"%ld %f %f %f \n",noeudsl[i].no,noeudsl[i].x,noeudsl[i].y,noeudsl[i].z);
-    i++;
- }
- fprintf (ft,"$ENDNOD\n");
-
-/* Ecriture des elements 
-########################*/
- fprintf (ft,"$ELM\n");
- fprintf (ft,"%ld\n",nbE);
- 
- i=0;
- while (i <nbE)
- {
-    if (elemsl[i].type == 2 ){
-       fprintf (ft,"%ld %ld %ld 1 3 %ld %ld %ld  \n",elemsl[i].n,elemsl[i].type,elemsl[i].grp,elemsl[i].n1,elemsl[i].n2,elemsl[i].n3);
-    }
-    else if (elemsl[i].type == 4 ){
-       fprintf (ft,"%ld %ld %ld %ld 4 %ld %ld %ld  %ld\n",elemsl[i].n,elemsl[i].type,elemsl[i].grp,elemsl[i].grp,elemsl[i].n1,elemsl[i].n2,elemsl[i].n3,elemsl[i].n4);
-    }
-    i++;
- }
- 
- fprintf (ft,"$ENDELM\n");
- 
- fclose(ft);
- printf ("Fichier Gmsh:%s ecrit !!\n",outdata);
-
- return 0;
-}
-/*#############################################################*/
-
-/*#############################################################*/
-int main(int argc,char *argv[])
-{
-
- datain();
-
- dataout();
-
- return 0;
-}