diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e53dbd37a4e61955c08d7ae4adfd2dce1b640e26..2450a74e37f3bf80653a2ac34be7ad91165a6480 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1237,44 +1237,6 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
   return true;
 }
 
-MElement *GModel::_createElementMSH(int num, int typeMSH, int physical,
-                                    int reg, int part, std::vector<MVertex*> &v,
-                                    std::map<int, std::vector<MElement*> > elements[10],
-                                    std::map<int, std::map<int, std::string> > physicals[4])
-{
-  if(CTX::instance()->mesh.switchElementTags) {
-    int tmp = reg;
-    reg = physical;
-    physical = tmp;
-  }
-
-  MElementFactory factory;
-  MElement *e = factory.create(typeMSH, v, num, part);
-  if(!e){
-    Msg::Error("Unknown type of element %d", typeMSH);
-    return 0;
-  }
-
-  switch(e->getType()){
-  case TYPE_PNT : elements[0][reg].push_back(e); break;
-  case TYPE_LIN : elements[1][reg].push_back(e); break;
-  case TYPE_TRI : elements[2][reg].push_back(e); break;
-  case TYPE_QUA : elements[3][reg].push_back(e); break;
-  case TYPE_TET : elements[4][reg].push_back(e); break;
-  case TYPE_HEX : elements[5][reg].push_back(e); break;
-  case TYPE_PRI : elements[6][reg].push_back(e); break;
-  case TYPE_PYR : elements[7][reg].push_back(e); break;
-  default : Msg::Error("Wrong type of element"); return 0;
-  }
-
-  int dim = e->getDim();
-  if(physical && (!physicals[dim].count(reg) || !physicals[dim][reg].count(physical)))
-    physicals[dim][reg][physical] = "unnamed";
-
-  if(part) getMeshPartitions().insert(part);
-  return e;
-}
-
 GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
                              std::vector<int> &elementNum,
                              std::vector<std::vector<int> > &vertexIndices,
@@ -1283,24 +1245,34 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
                              std::vector<int> &elementary,
                              std::vector<int> &partition)
 {
-  GModel *gm = new GModel();
-  std::map<int, std::vector<MElement*> > elements[10];
-  std::map<int, std::map<int, std::string> > physicals[4];
-  std::vector<MVertex*> vertexVector;
-
   int numVertices = (int)vertexMap.size();
   int numElement = (int)elementNum.size();
 
-  if(numElement != (int)vertexIndices.size())
+  if(numElement != (int)vertexIndices.size()){
     Msg::Error("Dimension in vertices numbers");
-  if(numElement != (int)elementType.size())
+    return 0;
+  }
+  if(numElement != (int)elementType.size()){
     Msg::Error("Dimension in elementType numbers");
-  if(numElement != (int)physical.size())
+    return 0;
+  }
+  if(numElement != (int)physical.size()){
     Msg::Error("Dimension in physical numbers");
-  if(numElement != (int)elementary.size())
+    return 0;
+  }
+  if(numElement != (int)elementary.size()){
     Msg::Error("Dimension in elementary numbers");
-  if(numElement != (int)partition.size())
+    return 0;
+  }
+  if(numElement != (int)partition.size()){
     Msg::Error("Dimension in partition numbers");
+    return 0;
+  }
+
+  GModel *gm = new GModel();
+  std::map<int, std::vector<MElement*> > elements[10];
+  std::map<int, std::map<int, std::string> > physicals[4];
+  std::vector<MVertex*> vertexVector;
 
   std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
   std::map<int, MVertex*>::const_iterator end = vertexMap.end();
@@ -1317,8 +1289,8 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
   if(minVertex == std::numeric_limits<int>::max())
     Msg::Error("Could not determine the min index of vertices");
 
-  // If the vertex numbering is dense, transfer the map into a
-  // vector to speed up element creation
+  // if the vertex numbering is dense, transfer the map into a vector to speed
+  // up element creation
   if((minVertex == 1 && maxVertex == numVertices) ||
      (minVertex == 0 && maxVertex == numVertices - 1)){
     Msg::Info("Vertex numbering is dense");
@@ -1355,8 +1327,31 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
       }
     }
 
-    gm->_createElementMSH(num, elementType[i], physical[i], elementary[i],
-                          partition[i], vertices, elements, physicals);
+    MElementFactory f;
+    MElement *e = f.create(elementType[i], vertices, num, partition[i]);
+    if(!e){
+      Msg::Error("Unknown type of element %d", elementType[i]);
+      delete gm;
+      return 0;
+    }
+    switch(e->getType()){
+    case TYPE_PNT : elements[0][elementary[i]].push_back(e); break;
+    case TYPE_LIN : elements[1][elementary[i]].push_back(e); break;
+    case TYPE_TRI : elements[2][elementary[i]].push_back(e); break;
+    case TYPE_QUA : elements[3][elementary[i]].push_back(e); break;
+    case TYPE_TET : elements[4][elementary[i]].push_back(e); break;
+    case TYPE_HEX : elements[5][elementary[i]].push_back(e); break;
+    case TYPE_PRI : elements[6][elementary[i]].push_back(e); break;
+    case TYPE_PYR : elements[7][elementary[i]].push_back(e); break;
+    case TYPE_POLYG: elements[8][elementary[i]].push_back(e); break;
+    case TYPE_POLYH: elements[9][elementary[i]].push_back(e); break;
+    default : Msg::Error("Wrong type of element"); delete gm; return 0;
+    }
+    int dim = e->getDim();
+    if(physical[i] && (!physicals[dim].count(elementary[i]) ||
+                       !physicals[dim][elementary[i]].count(physical[i])))
+      physicals[dim][elementary[i]][physical[i]] = "unnamed";
+    if(partition[i]) gm->getMeshPartitions().insert(partition[i]);
   }
 
   // store the elements in their associated elementary entity. If the
diff --git a/Geo/GModel.h b/Geo/GModel.h
index ea9b32a5e904fc0641b933b2cdf41b763d69b39c..f8f97c37be9b78b502eb4354b15f95ed2b9f83bb 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -123,13 +123,6 @@ class GModel
   void _storePhysicalTagsInEntities(int dim,
                                     std::map<int, std::map<int, std::string> > &map);
 
-  // utility function to create a mesh element in the I/O routines (don't use in
-  // new code; it will be removed eventually)
-  MElement *_createElementMSH(int num, int typeMSH, int physical,
-                              int reg, int part, std::vector<MVertex*> &v,
-                              std::map<int, std::vector<MElement*> > elements[10],
-                              std::map<int, std::map<int, std::string> > physicals[4]);
-
   // entity that is currently being meshed (used for error reporting)
   GEntity *_currentMeshEntity;
 
diff --git a/Geo/GModelIO_DIFF.cpp b/Geo/GModelIO_DIFF.cpp
index 0e5310f061a506a835233e0b5e31642d38e42aad..556e589be561496314f14acbe98eec7238addd33 100644
--- a/Geo/GModelIO_DIFF.cpp
+++ b/Geo/GModelIO_DIFF.cpp
@@ -294,9 +294,31 @@ int GModel::readDIFF(const std::string &name)
         if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices))
           return 0;
       }
-      _createElementMSH(num, type, physical, elementary[i-1][1], partition,
-                        vertices, elements, physicals);
-      // trouble if elementary[i-1][0]>1 nodal post-processing needed ?
+
+      MElementFactory f;
+      MElement *e = f.create(type, vertices, num, partition);
+      if(!e){
+        Msg::Error("Unknown type of element %d", type);
+        return 0;
+      }
+      int reg = elementary[i-1][1];
+      switch(e->getType()){
+      case TYPE_PNT : elements[0][reg].push_back(e); break;
+      case TYPE_LIN : elements[1][reg].push_back(e); break;
+      case TYPE_TRI : elements[2][reg].push_back(e); break;
+      case TYPE_QUA : elements[3][reg].push_back(e); break;
+      case TYPE_TET : elements[4][reg].push_back(e); break;
+      case TYPE_HEX : elements[5][reg].push_back(e); break;
+      case TYPE_PRI : elements[6][reg].push_back(e); break;
+      case TYPE_PYR : elements[7][reg].push_back(e); break;
+      default : Msg::Error("Wrong type of element"); return 0;
+      }
+      int dim = e->getDim();
+      if(physical && (!physicals[dim].count(reg) ||
+                      !physicals[dim][reg].count(physical)))
+        physicals[dim][reg][physical] = "unnamed";
+      if(partition) getMeshPartitions().insert(partition);
+
       if(numElements > 100000)
         Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
     }