diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index 2a067ad69c0dd67eec0f3ea6f7eca0d47bf92f83..cce4eb6d8b096dde3795fa645256d0f89bf94812 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_MED.cpp,v 1.21 2008-03-30 17:45:12 geuzaine Exp $
+// $Id: GModelIO_MED.cpp,v 1.22 2008-03-30 20:45:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -59,6 +59,67 @@ static int getElementTypeForMED(int msh, med_geometrie_element &med)
   }
 }
 
+int med2msh(med_geometrie_element med, int k)
+{
+  switch(med) {
+  case MED_SEG2: return k;
+  case MED_TRIA3: {
+    static const int map[3] = {0, 2, 1}; 
+    return map[k]; 
+  }
+  case MED_QUAD4: { 
+    static const int map[4] = {0, 3, 2, 1};
+    return map[k];
+  }
+  case MED_TETRA4: {
+    static const int map[4] = {0, 2, 1, 3};
+    return map[k];
+  }
+  case MED_HEXA8: {
+    static const int map[8] = {0, 3, 2, 1, 4, 7, 6, 5};
+    return map[k];
+  }
+  case MED_PENTA6: {
+    static const int map[6] = {0, 2, 1, 3, 5, 4};
+    return map[k];
+  }
+  case MED_PYRA5: {
+    static const int map[5] = {0, 3, 2, 1, 4};
+    return map[k];
+  }
+  case MED_SEG3: return k;
+  case MED_TRIA6: {
+    static const int map[6] = {0, 2, 1, 5, 4, 3};
+    return map[k];
+  }
+  case MED_TETRA10: {
+    static const int map[10] = {0, 2, 1, 3, 6, 5, 4, 7, 8, 9};
+    return map[k];
+  }
+  case MED_POINT1: return k;
+  case MED_QUAD8: {
+    static const int map[8] = {0, 3, 2, 1, 7, 6, 5, 4};
+    return map[k];
+  }
+  case MED_HEXA20: {
+    static const int map[20] = {0, 3, 2, 1, 4, 7, 6, 5, 11, 8, 16,
+				10, 19, 9, 18, 17, 15, 12, 14, 13};
+    return map[k];
+  }
+  case MED_PENTA15: {
+    static const int map[15] = {0, 2, 1, 3, 5, 4, 8, 6, 12, 7, 14, 13, 11, 9, 10};
+    return map[k];
+  }
+  case MED_PYRA13: {
+    static const int map[13] = {0, 3, 2, 1, 4, 8, 5, 9, 7, 12, 6, 11, 10};
+    return map[k];
+  }
+  default:
+    Msg(GERROR, "Unknown MED element type");
+    return k;
+  }
+}
+
 int GModel::readMED(const std::string &name)
 {
   med_idt fid = MEDouvrir((char*)name.c_str(), MED_LECTURE);
@@ -180,7 +241,7 @@ int GModel::readMED(const std::string &name, int meshIndex)
       for(int j = 0; j < numEle; j++){    
 	std::vector<MVertex*> v(numNodPerEle);
 	for(int k = 0; k < numNodPerEle; k++)
-	  v[k] = verts[conn[numNodPerEle * j + k] - 1];
+	  v[k] = verts[conn[numNodPerEle * j + med2msh(type, k)] - 1];
 	MElement *e = factory.create(mshType, v, eleTags.empty() ? 0 : eleTags[j]);
 	if(e) elements[-fam[j]].push_back(e);
       }
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index a16c539e20c56bb2243c7161685ea4fb93cc2fa8..0402c32039e5f5033309d30b6f3715aefc6b4443 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.48 2008-03-29 21:36:29 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.49 2008-03-30 20:45:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -1810,7 +1810,10 @@ int GModel::readBDF(const std::string &name)
       else if(!strncmp(buffer, "CTETRA", 6)){
         if(readElementBDF(fp, buffer, 6, -4, num, region, vertices, vertexMap)){
           if(vertices.size() == 10)
-            elements[3][region].push_back(new MTetrahedron10(vertices, num));
+            elements[3][region].push_back
+	      (new MTetrahedron10(vertices[0], vertices[1], vertices[2], vertices[3], 
+				  vertices[4], vertices[5], vertices[6], vertices[7], 
+				  vertices[9], vertices[8], num));
           else
             elements[3][region].push_back(new MTetrahedron(vertices, num));
         }
@@ -1818,7 +1821,13 @@ int GModel::readBDF(const std::string &name)
       else if(!strncmp(buffer, "CHEXA", 5)){
         if(readElementBDF(fp, buffer, 5, -8, num, region, vertices, vertexMap)){
           if(vertices.size() == 20)
-            elements[4][region].push_back(new MHexahedron20(vertices, num));
+            elements[4][region].push_back
+	      (new MHexahedron20(vertices[0], vertices[1], vertices[2], vertices[3], 
+				 vertices[4], vertices[5], vertices[6], vertices[7], 
+				 vertices[8], vertices[11], vertices[12], vertices[9], 
+				 vertices[13], vertices[10], vertices[14], vertices[15], 
+				 vertices[16], vertices[19], vertices[17], vertices[18], 
+				 num));
           else
             elements[4][region].push_back(new MHexahedron(vertices, num));
         }
@@ -1826,7 +1835,11 @@ int GModel::readBDF(const std::string &name)
       else if(!strncmp(buffer, "CPENTA", 6)){
         if(readElementBDF(fp, buffer, 6, -6, num, region, vertices, vertexMap)){
           if(vertices.size() == 15)
-            elements[5][region].push_back(new MPrism15(vertices, num));
+            elements[5][region].push_back
+	      (new MPrism15(vertices[0], vertices[1], vertices[2], vertices[3],
+			    vertices[4], vertices[5], vertices[6], vertices[8],
+			    vertices[9], vertices[7], vertices[10], vertices[11],
+			    vertices[12], vertices[14], vertices[13], num));
           else
             elements[5][region].push_back(new MPrism(vertices, num));
         }
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 38ac7a2fd5e14018c876708f973ac3f34e6ece18..56b7fe3f2e7c994eb6143b8e561d6df91cbf0367 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.63 2008-03-29 21:36:29 geuzaine Exp $
+// $Id: MElement.cpp,v 1.64 2008-03-30 20:45:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -583,7 +583,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   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)->getIndex());
+      fprintf(fp, ",%d", getVertexBDF(i)->getIndex());
       if(i != n - 1 && !((i + 3) % 8)){
         fprintf(fp, ",+%s%d\n+%s%d", cont[ncont], _num, cont[ncont], _num);
         ncont++;
@@ -596,7 +596,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
   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)->getIndex());
+      fprintf(fp, "%-8d", getVertexBDF(i)->getIndex());
       if(i != n - 1 && !((i + 3) % 8)){
         fprintf(fp, "+%s%-6d\n+%s%-6d", cont[ncont], _num, cont[ncont], _num);
         ncont++;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 90b654f9c3fa0789a2ce5de4ec9ba40dfb0b2878..1e3130306d3d486fa88d723654b11e6d219c05cf 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -1038,6 +1038,11 @@ class MHexahedron : public MElement {
   virtual int getDim(){ return 3; }
   virtual int getNumVertices(){ return 8; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[8] = {0, 3, 2, 1, 4, 7, 6, 5};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdges(){ return 12; }
   virtual MEdge getEdge(int num)
   {
@@ -1196,6 +1201,12 @@ class MHexahedron20 : public MHexahedron {
                                 9, 10, 12, 14, 15, 16, 18, 19, 17};
     return getVertex(map[num]); 
   }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[20] = {0, 3, 2, 1, 4, 7, 6, 5, 9, 13, 11, 
+				8, 17, 19, 18, 16, 10, 15, 14, 12};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 12; }
   virtual int getNumEdgesRep(){ return 24; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1364,6 +1375,11 @@ class MPrism : public MElement {
   virtual int getDim(){ return 3; }
   virtual int getNumVertices(){ return 6; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[6] = {0, 2, 1, 3, 5, 4};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdges(){ return 9; }
   virtual MEdge getEdge(int num)
   {
@@ -1512,6 +1528,11 @@ class MPrism15 : public MPrism {
     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 MVertex *getVertexMED(int num)
+  {
+    static const int map[15] = {0, 2, 1, 3, 5, 4, 7, 9, 6, 13, 14, 12, 8, 11, 10};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 9; }
   virtual int getNumEdgesRep(){ return 18; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1655,6 +1676,11 @@ class MPyramid : public MElement {
   virtual int getDim(){ return 3; }
   virtual int getNumVertices(){ return 5; }
   virtual MVertex *getVertex(int num){ return _v[num]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[5] = {0, 3, 2, 1, 4};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdges(){ return 8; }
   virtual MEdge getEdge(int num)
   {
@@ -1789,6 +1815,11 @@ class MPyramid13 : public MPyramid {
   virtual int getPolynomialOrder(){ return 2; }
   virtual int getNumVertices(){ return 13; }
   virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
+  virtual MVertex *getVertexMED(int num)
+  {
+    static const int map[13] = {0, 3, 2, 1, 4, 6, 10, 8, 5, 7, 12, 11, 9};
+    return getVertex(map[num]); 
+  }
   virtual int getNumEdgeVertices(){ return 8; }
   virtual int getNumEdgesRep(){ return 16; }
   virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index d3b1e71d94e6b660d39e02d91a212d26dcc07007..39a673cdf393e783267d470086c132c9f5681d78 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModelIO.cpp,v 1.25 2008-03-30 17:45:12 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.26 2008-03-30 20:45:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -321,9 +321,9 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	  num = tags[profile[i] - 1];
 	}
 	double *d = _steps[step]->getData(num, true, mult);
-	// FIXME: for ElementNodeData, we need to reorder the data before
-	// storing it
-	// FIXME: what should we do with Gauss point data?
+	// FIXME: for ElementNodeData, we need to reorder the data
+	// before storing it, using med2msh()). Also: what should we
+	// do with Gauss point data?
 	for(int j = 0; j < numComp * mult; j++)
 	  d[j] = val[numComp * i + j];
 	double s = ComputeScalarRep(_steps[step]->getNumComponents(), d);
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 751d6b8720ddce7f27ffe41ee24711dc977ad45e..c8fb22f258d56c002580e059ff6fb80a5124d472 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,7 +1,8 @@
-$Id: VERSIONS,v 1.401 2008-03-21 17:09:06 geuzaine Exp $
+$Id: VERSIONS,v 1.402 2008-03-30 20:45:28 geuzaine Exp $
 
-2.1.2 (xx): new Fields interface; more work on new post-processing
-API; fixed various bugs.
+2.1.2 (xx): new Fields interface; new model-based post-processing
+backend; added MED I/O for mesh and post-processing; fixed BDF vertex
+ordering for 2nd order elements; fixed various bugs.
 
 2.1.1 (Mar 1, 2008): small bug fixes (second order meshes, combine
 views, divide and conquer crash, ...).