diff --git a/Geo/GModelIO_MESH.cpp b/Geo/GModelIO_MESH.cpp
index 6e1a70bb002793a133881449d91079e48fe4f618..deca5660c5b6c7a4b7bb4fd6cc9b548d4ee32bd1 100644
--- a/Geo/GModelIO_MESH.cpp
+++ b/Geo/GModelIO_MESH.cpp
@@ -11,6 +11,7 @@
 #include "MQuadrangle.h"
 #include "MTetrahedron.h"
 #include "MHexahedron.h"
+#include "Context.h"
 
 static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
                         std::vector<MVertex*> &vertices)
@@ -100,6 +101,21 @@ int GModel::readMESH(const std::string &name)
           elements[1][cl].push_back(new MTriangle(vertices));
         }
       }
+      else if(!strcmp(str, "TrianglesP2")){
+        if(!fgets(buffer, sizeof(buffer), fp)) break;
+        int nbe;
+        sscanf(buffer, "%d", &nbe);
+        Msg::Info("%d triangles", nbe);
+        for(int i = 0; i < nbe; i++) {
+          if(!fgets(buffer, sizeof(buffer), fp)) break;
+          int n[6], cl;
+          sscanf(buffer, "%d %d %d %d", &n[0], &n[1], &n[2], &n[3], &n[4], &n[5],&cl);
+          for(int j = 0; j < 6; j++) n[j]--;
+          std::vector<MVertex*> vertices;
+          if(!getVertices(6, n, vertexVector, vertices)) return 0;
+          elements[1][cl].push_back(new MTriangle6(vertices));
+        }
+      }
       else if(!strcmp(str, "Quadrilaterals")) {
         if(!fgets(buffer, sizeof(buffer), fp)) break;
         int nbe;
@@ -130,6 +146,21 @@ int GModel::readMESH(const std::string &name)
           elements[3][cl].push_back(new MTetrahedron(vertices));
         }
       }
+      else if(!strcmp(str, "TetrahedraP2")) {
+        if(!fgets(buffer, sizeof(buffer), fp)) break;
+        int nbe;
+        sscanf(buffer, "%d", &nbe);
+        Msg::Info("%d tetrahedra", nbe);
+        for(int i = 0; i < nbe; i++) {
+          if(!fgets(buffer, sizeof(buffer), fp)) break;
+          int n[10], cl;
+          sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3],&n[4], &n[5], &n[6], &n[7], &n[9], &n[8], &cl);
+          for(int j = 0; j < 10; j++) n[j]--;
+          std::vector<MVertex*> vertices;
+          if(!getVertices(10, n, vertexVector, vertices)) return 0;
+          elements[3][cl].push_back(new MTetrahedron10(vertices));
+        }
+      }
       else if(!strcmp(str, "Hexahedra")) {
         if(!fgets(buffer, sizeof(buffer), fp)) break;
         int nbe;
@@ -216,7 +247,11 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
     }
   }
   if(numTriangles){
-    fprintf(fp, " Triangles\n");
+    if(CTX::instance()->mesh.order == 2)
+      fprintf(fp, " TrianglesP2\n");
+    else
+      fprintf(fp, " Triangles\n");
+
     fprintf(fp, " %d\n", numTriangles);
     for(fiter it = firstFace(); it != lastFace(); ++it){
       int numPhys = (*it)->physicals.size();
@@ -240,7 +275,10 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
     }
   }
   if(numTetrahedra){
-    fprintf(fp, " Tetrahedra\n");
+    if(CTX::instance()->mesh.order == 2)
+      fprintf(fp, " TetrahedraP2\n");
+    else
+      fprintf(fp, " Tetrahedra\n");
     fprintf(fp, " %d\n", numTetrahedra);
     for(riter it = firstRegion(); it != lastRegion(); ++it){
       int numPhys = (*it)->physicals.size();