From 9dcabb1e055104809fc82f8eb0a5c22f419fd3b4 Mon Sep 17 00:00:00 2001
From: Stefen Guzik <guzik2@llnl.gov>
Date: Tue, 5 Aug 2008 07:28:40 +0000
Subject: [PATCH] Code simplification

---
 Geo/ElementTraits.h                | 176 -----------------
 Geo/GEdge.cpp                      |  10 +
 Geo/GEdge.h                        |  11 +-
 Geo/GEntity.h                      |   9 +-
 Geo/GFace.cpp                      |  16 ++
 Geo/GFace.h                        |  11 +-
 Geo/GModel.cpp                     |  15 ++
 Geo/GModel.h                       |  11 ++
 Geo/GModelIO_CGNS.cpp              |  35 +++-
 Geo/GRegion.cpp                    |  22 +++
 Geo/GRegion.h                      |  11 +-
 Geo/MElement.h                     |  38 ++++
 Geo/MZone.cpp                      | 145 +++++---------
 Geo/MZone.h                        |   1 -
 Geo/Makefile                       |   6 +-
 Mesh/Makefile                      |  14 +-
 Mesh/Partition.cpp                 | 307 ++++++++++-------------------
 contrib/Chaco/main/Gmsh_printf.cpp |   1 +
 18 files changed, 351 insertions(+), 488 deletions(-)
 delete mode 100644 Geo/ElementTraits.h

diff --git a/Geo/ElementTraits.h b/Geo/ElementTraits.h
deleted file mode 100644
index 7252ec3681..0000000000
--- a/Geo/ElementTraits.h
+++ /dev/null
@@ -1,176 +0,0 @@
-// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-//
-// ElementTraits.h - Copyright (C) 2008 S. Guzik, C. Geuzaine, J.-F. Remacle
-#ifndef _ELEMENTTRAITS_H_
-#define _ELEMENTTRAITS_H_
-
-#include "MElement.h"
-#include "GEdge.h"
-#include "GFace.h"
-#include "GRegion.h"
-
-
-/*******************************************************************************
- *
- * - The classes in this file define traits for elements
- *
- ******************************************************************************/
-
-//--Traits based on the dimension
-
-template <unsigned DIM> struct DimTr;
-template <> struct DimTr<1>
-{
-  typedef GEdge EntityT;
-  enum { numElemTypes = 1 };
-};
-template <> struct DimTr<2>
-{
-  typedef MEdge FaceT;
-  typedef GFace EntityT;
-  enum { numElemTypes = 2 };
-};
-template <> struct DimTr<3>
-{
-  typedef MFace FaceT;
-  typedef GRegion EntityT;
-  enum { numElemTypes = 4 };
-};
-
-//--This traits class describes the number of dimension-based 'FaceT' in an
-//--primary element type
-
-template <unsigned DIM, typename Elem> struct ElemFaceTr;
-template <> struct ElemFaceTr<2, MLine>        { enum { numFaceT =  1 }; };
-template <> struct ElemFaceTr<2, MTriangle>    { enum { numFaceT =  3 }; };
-template <> struct ElemFaceTr<2, MQuadrangle>  { enum { numFaceT =  4 }; };
-template <> struct ElemFaceTr<3, MTriangle>    { enum { numFaceT =  1 }; };
-template <> struct ElemFaceTr<3, MQuadrangle>  { enum { numFaceT =  1 }; };
-template <> struct ElemFaceTr<3, MTetrahedron> { enum { numFaceT =  4 }; };
-template <> struct ElemFaceTr<3, MHexahedron>  { enum { numFaceT =  6 }; };
-template <> struct ElemFaceTr<3, MPrism>       { enum { numFaceT =  5 }; };
-template <> struct ElemFaceTr<3, MPyramid>     { enum { numFaceT =  5 }; };
-
-//--This traits class gives iterator types and begin and end iterators for
-//--element type number N in an entity of dimension DIM
-
-template <unsigned DIM, int N> struct DimElemTr;
-template <> struct DimElemTr<1, 1> {
-  typedef MLine Elem;
-  typedef std::vector<MLine*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GEdge *const gEdge)
-  {
-    return gEdge->lines.begin();
-  }
-  static ConstElementIterator end(const GEdge *const gEdge)
-  {
-    return gEdge->lines.end();
-  }
-};
-template <> struct DimElemTr<2, 1> {
-  typedef MQuadrangle Elem;
-  typedef std::vector<MQuadrangle*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GFace *const gFace)
-  {
-    return gFace->quadrangles.begin();
-  }
-  static ConstElementIterator end(const GFace *const gFace)
-  {
-    return gFace->quadrangles.end();
-  }
-};
-template <> struct DimElemTr<2, 2> {
-  typedef MTriangle Elem;
-  typedef std::vector<MTriangle*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GFace *const gFace)
-  {
-    return gFace->triangles.begin();
-  }
-  static ConstElementIterator end(const GFace *const gFace)
-  {
-    return gFace->triangles.end();
-  }
-};
-template <> struct DimElemTr<3, 1> {
-  typedef MPyramid Elem;
-  typedef std::vector<MPyramid*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GRegion *const gRegion)
-  {
-    return gRegion->pyramids.begin();
-  }
-  static ConstElementIterator end(const GRegion *const gRegion)
-  {
-    return gRegion->pyramids.end();
-  }
-};
-template <> struct DimElemTr<3, 2> {
-  typedef MPrism Elem;
-  typedef std::vector<MPrism*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GRegion *const gRegion)
-  {
-    return gRegion->prisms.begin();
-  }
-  static ConstElementIterator end(const GRegion *const gRegion)
-  {
-    return gRegion->prisms.end();
-  }
-};
-template <> struct DimElemTr<3, 3> {
-  typedef MHexahedron Elem;
-  typedef std::vector<MHexahedron*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GRegion *const gRegion)
-  {
-    return gRegion->hexahedra.begin();
-  }
-  static ConstElementIterator end(const GRegion *const gRegion)
-  {
-    return gRegion->hexahedra.end();
-  }
-};
-template <> struct DimElemTr<3, 4> {
-  typedef MTetrahedron Elem;
-  typedef std::vector<MTetrahedron*>::const_iterator ConstElementIterator;
-  static ConstElementIterator begin(const GRegion *const gRegion)
-  {
-    return gRegion->tetrahedra.begin();
-  }
-  static ConstElementIterator end(const GRegion *const gRegion)
-  {
-    return gRegion->tetrahedra.end();
-  }
-};
-
-//--Traits/policy class based on dimension of a face.
-
-template <typename FaceT> struct FaceTr;
-template <> struct FaceTr<MEdge>
-{
-  template <typename Elem>
-  static MEdge getFace(Elem *const element, const int iFace)
-  {
-    return element->getEdge(iFace);  // Primary element
-  }
-  static void getAllFaceVertices(MElement *const element, const int iFace,
-                                 std::vector<MVertex*> &v)
-  {
-    element->getEdgeVertices(iFace, v);  // Any element
-  }
-};
-template <> struct FaceTr<MFace>
-{
-  template <typename Elem>
-  static MFace getFace(Elem *const element, const int iFace)
-  {
-    return element->getFace(iFace);  // Primary element
-  }
-  static void getAllFaceVertices(MElement *const element, const int iFace,
-                                 std::vector<MVertex*> &v)
-  {
-    element->getFaceVertices(iFace, v);  // Any element
-  }
-};
-
-#endif
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index c83a83580a..ad92be02e4 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -42,6 +42,16 @@ unsigned int GEdge::getNumMeshElements()
   return lines.size();
 }
 
+void GEdge::getNumMeshElements(unsigned *const c) const
+{
+  c[0] += lines.size();
+}
+
+MElement *const *GEdge::getStartElementType(int type) const
+{
+  return reinterpret_cast<MElement *const *>(&lines[0]);
+}
+
 MElement *GEdge::getMeshElement(unsigned int index)
 { 
   if(index < lines.size())
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index ee7e6b7302..4c9a897805 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -107,8 +107,17 @@ class GEdge : public GEntity {
   // true if start == end and no more than 2 segments
   bool isMeshDegenerated() const{ return (v0 == v1 && mesh_vertices.size() < 2); }
 
-  // get number of elements in the mesh and get element by index
+  // number of types of elements
+  int getNumElementTypes() const { return 1; }
+
+  // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  void getNumMeshElements(unsigned *const c) const;
+
+  // get the start of the array of a type of element
+  MElement *const *getStartElementType(int type) const;
+
+  // get the element at the given index
   MElement *getMeshElement(unsigned int index);
 
   // reset the mesh attributes to default values
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 3596aa28a5..5476df2cc7 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -223,8 +223,15 @@ class GEntity {
   // reset the mesh attributes to default values
   virtual void resetMeshAttributes() { return; }
 
-  // get the number of mesh elements in the entity
+  // number of types of elements
+  virtual int getNumElementTypes() const { return 0; }
+
+  // get the number of mesh elements (total and by type) in the entity
   virtual unsigned int getNumMeshElements() { return 0; }
+  virtual void getNumMeshElements(unsigned *const c) const { };
+
+  // get the start of the array of a type of element
+  virtual MElement *const *getStartElementType(int type) const { return 0; }
 
   // get the element at the given index
   virtual MElement *getMeshElement(unsigned int index) { return 0; }
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index f20e23be58..4b155a4ac4 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -64,6 +64,22 @@ unsigned int GFace::getNumMeshElements()
   return triangles.size() + quadrangles.size(); 
 }
 
+void GFace::getNumMeshElements(unsigned *const c) const
+{
+  c[0] += triangles.size();
+  c[1] += quadrangles.size();
+}
+
+MElement *const *GFace::getStartElementType(int type) const
+{
+  switch(type) {
+  case 0:
+    return reinterpret_cast<MElement *const *>(&triangles[0]);
+  case 1:
+    return reinterpret_cast<MElement *const *>(&quadrangles[0]);
+  }
+}
+
 MElement *GFace::getMeshElement(unsigned int index)
 { 
   if(index < triangles.size())
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 37c10e92ae..f20f3d9f09 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -154,8 +154,17 @@ class GFace : public GEntity
                         double &x, double &y, double &z) const;
   void getMeanPlaneData(double plan[3][3]) const;
 
-  // get number of elements in the mesh and get element by index
+  // number of types of elements
+  int getNumElementTypes() const { return 2; }
+
+  // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  void getNumMeshElements(unsigned *const c) const;
+
+  // get the start of the array of a type of element
+  MElement *const *getStartElementType(int type) const;
+
+  // get the element at the given index
   MElement *getMeshElement(unsigned int index);
 
   // reset the mesh attributes to default values
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 1a7daf14e7..a18bdd2f0e 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -433,6 +433,21 @@ int GModel::getNumMeshElements()
   return n;
 }
 
+int GModel::getNumMeshElements(unsigned c[4])
+{
+  c[0] = 0; c[1] = 0; c[2] = 0; c[3] = 0;
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    (*it)->getNumMeshElements(c);
+  if(c[0] + c[1] + c[2] + c[3]) return 3;
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+    (*it)->getNumMeshElements(c);
+  if(c[0] + c[1]) return 2;
+  for(eiter it = firstEdge(); it != lastEdge(); ++it)
+    (*it)->getNumMeshElements(c);
+  if(c[0]) return 1;
+  return 0;
+}
+
 static void MElementBB(void *a, double *min, double *max)
 {
   MElement *e = (MElement*)a;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 6a20badac8..f51b3acaf8 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -77,6 +77,7 @@ class GModel
   std::set<GVertex*, GEntityLessThan> vertices;
   std::set<int> meshPartitions;
   std::map<int, std::string> physicalNames, elementaryNames;
+  int partitionSize[2];
 
  public:
   GModel(std::string name="");
@@ -192,6 +193,10 @@ class GModel
   // Returns the total number of elements in the mesh
   int getNumMeshElements();
 
+  // Get the number of each type of element in the mesh at the largest
+  // dimension and return the dimension
+  int getNumMeshElements(unsigned c[4]);
+
   // Access a mesh element by coordinates
   MElement *getMeshElementByCoord(SPoint3 &p);
 
@@ -226,6 +231,12 @@ class GModel
   // Deletes all the partitions
   void deleteMeshPartitions();
 
+  // Store/recall min and max partitions size
+  void setMinPartitionSize(const int pSize) { partitionSize[0] = pSize; }
+  void setMaxPartitionSize(const int pSize) { partitionSize[1] = pSize; }
+  int getMinPartitionSize() const { return partitionSize[0]; }
+  int getMaxPartitionSize() const { return partitionSize[1]; }
+
   // Performs various coherence tests on the mesh
   void checkMeshCoherence();
 
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 8903a96862..7e589938e6 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -216,7 +216,7 @@ int GModel::readCGNS(const std::string &name)
  *
  ******************************************************************************/
 
-int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
+int GModel::writeCGNS(const std::string &name, int zoneDefinition,
                       const CGNSOptions &options, double scalingFactor)
 {
 
@@ -230,13 +230,44 @@ int GModel::writeCGNS(const std::string &name, const int zoneDefinition,
   PhysGroupMap groups[4];               // vector of entities that belong to
                                         // each physical group (in each
                                         // dimension)
+  int numZone;
+  std::vector<std::vector<MElement*> > zoneElements;
+  switch(zoneDefinition) {
+  case 1:  // By partition
+    {
+      numZone = meshPartitions.size();
+      zoneElements.resize(numZone);
+      for(int i = 0; i != numZone; ++i)
+        zoneElements[i].reserve(getMaxPartitionSize());
+//       typedef GRegion Ent;
+//       Ent *ent;
+//       unsigned numElem[4];
+//       numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+//       ent->getNumMeshElements(numElem);
+//       int nType = ent->getNumTypeElements();
+//       for(int iType = 0; iType != nType; ++iType) {
+//         MElement *const start = ent->getMeshElement
+//           ((iType == 0) ? 0 : numElem[iType-1]);
+//         // Then can loop
+//         const int nElem = numElem[iType];
+//         for(int iElem = 0; iElem != nElem; ++iElem) {
+//         }
+//       }
+
+    }
+    break;
+  case 2:  // By physical
+     break;
+  }
 
 //--Get a list of groups in each dimension and each of the entities in that
 //--group.  The groups will define the zones.  If no 2D or 3D groups exist, just
 //--store all mesh elements in one zone.
 
   getPhysicalGroups(groups);
-  if(groups[face].size() + groups[region].size() == 0) {
+  if(groups[face].size() + groups[region].size() == 0) zoneDefinition = 0;
+  {
+
     // If there are no 2D or 3D physical groups, save all elements in one zone.
     // Pretend that there is one physical group ecompassing all the elements.
     for(riter it = firstRegion(); it != lastRegion(); ++it)
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index ced3ab2415..b4f55b540d 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -48,6 +48,28 @@ unsigned int GRegion::getNumMeshElements()
   return tetrahedra.size() + hexahedra.size() + prisms.size() + pyramids.size();
 }
 
+void GRegion::getNumMeshElements(unsigned *const c) const
+{
+  c[0] += tetrahedra.size();
+  c[1] += hexahedra.size();
+  c[3] += prisms.size();
+  c[4] += pyramids.size();
+}
+
+MElement *const *GRegion::getStartElementType(int type) const
+{
+  switch(type) {
+  case 0:
+    return reinterpret_cast<MElement *const *>(&tetrahedra[0]);
+  case 1:
+    return reinterpret_cast<MElement *const *>(&hexahedra[0]);
+  case 2:
+    return reinterpret_cast<MElement *const *>(&prisms[0]);
+  case 3:
+    return reinterpret_cast<MElement *const *>(&pyramids[0]);
+  }
+}
+
 MElement *GRegion::getMeshElement(unsigned int index)
 { 
   if(index < tetrahedra.size())
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 0368336169..ccb61be903 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -53,8 +53,17 @@ class GRegion : public GEntity {
   // return a type-specific additional information string
   virtual std::string getAdditionalInfoString();
 
-  // get number of elements in the mesh and get element by index
+  // number of types of elements
+  int getNumElementTypes() const { return 4; }
+
+  // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  void getNumMeshElements(unsigned *const c) const;
+
+  // get the start of the array of a type of element
+  MElement *const *getStartElementType(int type) const;
+
+  // get the element at the given index
   MElement *getMeshElement(unsigned int index);
 
   // reset the mesh attributes to default values
diff --git a/Geo/MElement.h b/Geo/MElement.h
index e8388230fb..2012099606 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -2850,4 +2850,42 @@ struct compareMTriangleLexicographic
   }
 };
 
+// Traits of various elements based on the dimension.  These generally define
+// the faces of 2-D elements as MEdge and 3-D elements as MFace.
+
+template <unsigned DIM> struct DimTr;
+template <> struct DimTr<2>
+{
+  typedef MEdge FaceT;
+  static int getNumFace(MElement *const element)
+  {
+    return element->getNumEdges();
+  }
+  static MEdge getFace(MElement *const element, const int iFace)
+  {
+    return element->getEdge(iFace);
+  }
+  static void getAllFaceVertices(MElement *const element, const int iFace,
+                                 std::vector<MVertex*> &v)
+  {
+    element->getEdgeVertices(iFace, v);
+  }
+};
+template <> struct DimTr<3>
+{
+  typedef MFace FaceT;
+  static int getNumFace(MElement *const element)
+  {
+    return element->getNumFaces();
+  }
+  static MFace getFace(MElement *const element, const int iFace)
+  {
+    return element->getFace(iFace);
+  }
+  static void getAllFaceVertices(MElement *const element, const int iFace,
+                                 std::vector<MVertex*> &v)
+  {
+    element->getFaceVertices(iFace, v);
+  }
+};
 #endif
diff --git a/Geo/MZone.cpp b/Geo/MZone.cpp
index 76f22e6901..6391e70b76 100644
--- a/Geo/MZone.cpp
+++ b/Geo/MZone.cpp
@@ -12,7 +12,7 @@
  * Forward declarations
  *============================================================================*/
 
-template <unsigned DIM, int N = DimTr<DIM>::numElemTypes>
+template <unsigned DIM>
 struct ParseEntity;
 
 
@@ -53,12 +53,10 @@ void MZone<DIM>::add_elements_in_entities
 (EntIter begin, EntIter end, const int partition)
 {
 
-  typedef typename DimTr<DIM>::EntityT Ent;
-
   // Find the neighbours of each vertex, edge, and face
   for(EntIter itEnt = begin; itEnt != end; ++itEnt) {
-    ParseEntity<DIM>::eval(static_cast<Ent*>(*itEnt), boFaceMap, elemVec,
-                           vertMap, zoneElemConn, partition);
+    ParseEntity<DIM>::eval(*itEnt, boFaceMap, elemVec, vertMap, zoneElemConn,
+                           partition);
   }
 
 }
@@ -96,11 +94,9 @@ void MZone<DIM>::add_elements_in_entity
 (EntPtr entity, const int partition)
 {
 
-  typedef typename DimTr<DIM>::EntityT Ent;
-
   // Find the neighbours of each vertex, edge, and face
-  ParseEntity<DIM>::eval(static_cast<Ent*>(entity), boFaceMap, elemVec,
-                         vertMap, zoneElemConn, partition);
+  ParseEntity<DIM>::eval(entity, boFaceMap, elemVec, vertMap, zoneElemConn,
+                         partition);
 
 }
 
@@ -136,7 +132,7 @@ int MZone<DIM>::zoneData()
   for(typename BoFaceMap::iterator fMapIt = boFaceMap.begin();
       fMapIt != boFaceMap.end(); ++fMapIt) {
     // Get all the vertices on the face
-    FaceTr<FaceT>::getAllFaceVertices
+    DimTr<DIM>::getAllFaceVertices
       (elemVec[fMapIt->second.parentElementIndex].element,
        fMapIt->second.parentFace, faceVertices);
     const int nVert = faceVertices.size();
@@ -247,108 +243,61 @@ int MZone<DIM>::zoneData()
 
 /*******************************************************************************
  *
- * Template meta-programming classes
+ * Struct ParseEntity
+ *
+ * Purpose
+ * =======
+ *
+ *   Iterates over the elements in an entity and adds them to the MZone.
  *
  ******************************************************************************/
 
-/*--------------------------------------------------------------------*
- * This class uses traits to determine the element types in an entity
- * and iterate over them.  A template meta-programming loop is used to
- * iterate over the types of elements.  So within the class, we are
- * simply examining one type of element of a specific dimension in one
- * specific entity.
- *--------------------------------------------------------------------*/
-
-//--Entry point
-
-template <unsigned DIM, int N>
+template <unsigned DIM>
 struct ParseEntity
 {
-  typedef typename DimTr<DIM>::EntityT Ent;  // The type of entity
   typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
   typedef typename LFaceTr<FaceT>::BoFaceMap BoFaceMap; // The corresponding map
-  typedef typename DimElemTr<DIM, N>::Elem Elem;  // The type of primary element
-  typedef typename DimElemTr<DIM, N>::ConstElementIterator ConstElementIterator;
-  static void eval(const Ent *const entity,
-                   BoFaceMap &boFaceMap,
-                   ElementVec &elemVec,
-                   VertexMap &vertMap,
-                   ElementConnectivity *zoneElemConn,
-                   const int partition)
-  {
-    ConstElementIterator elemEnd = DimElemTr<DIM, N>::end(entity);
-    for(ConstElementIterator itElem = DimElemTr<DIM, N>::begin(entity);
-        itElem != elemEnd; ++itElem) {
-      if(partition < 0 || (*itElem)->getPartition() == partition) {
-        // Unique list of elements
-        const unsigned eVecIndex = elemVec.size();
-        elemVec.push_back(ElemData(*itElem));
-        ++zoneElemConn[(*itElem)->getTypeForMSH() - 1].numElem;
-        // Unique list of vertices
-        const int nVert = (*itElem)->getNumVertices();
-        for(int iVert = 0; iVert != nVert; ++iVert)
-          vertMap[(*itElem)->getVertex(iVert)] = 0;  // Unlabelled
-        // Maintain list of (base element) faces with only one bounding element.
-        for(int iFace = 0; iFace != ElemFaceTr<DIM, Elem>::numFaceT;
-            ++iFace) {
-          FaceT face = FaceTr<FaceT>::template getFace<Elem>(*itElem, iFace);
-          std::pair<typename BoFaceMap::iterator, bool> insBoFaceMap =
-            boFaceMap.insert(std::pair<FaceT, FaceData>
-                             (face, FaceData(*itElem, iFace, eVecIndex)));
-          if(!insBoFaceMap.second) {
-            // The face already exists and is therefore bounded on both sides
-            // by elements.  Not a boundary face so delete.
-            boFaceMap.erase(insBoFaceMap.first);
-          }
-        }
-      }
-    }
-    // Next element type in the entity
-    ParseEntity<DIM, N-1>::eval(entity, boFaceMap, elemVec, vertMap,
-                                zoneElemConn, partition);
-  }
-};
-
-//--Terminate loop when N = 1
 
-template <unsigned DIM>
-struct ParseEntity<DIM, 1>
-{
-  typedef typename DimTr<DIM>::EntityT Ent;  // The type of entity
-  typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
-  typedef typename LFaceTr<FaceT>::BoFaceMap BoFaceMap; // The corresponding map
-  typedef typename DimElemTr<DIM, 1>::Elem Elem;  // The type of primary element
-  typedef typename DimElemTr<DIM, 1>::ConstElementIterator ConstElementIterator;
-  static void eval(const Ent *const entity,
+  static void eval(const GEntity *const entity,
                    BoFaceMap &boFaceMap,
                    ElementVec &elemVec,
                    VertexMap &vertMap,
                    ElementConnectivity *zoneElemConn,
                    const int partition)
   {
-    ConstElementIterator elemEnd = DimElemTr<DIM, 1>::end(entity);
-    for(ConstElementIterator itElem = DimElemTr<DIM, 1>::begin(entity);
-        itElem != elemEnd; ++itElem) {
-      if(partition < 0 || (*itElem)->getPartition() == partition) {
-        // Unique list of elements
-        const unsigned eVecIndex = elemVec.size();
-        elemVec.push_back(ElemData(*itElem));
-        ++zoneElemConn[(*itElem)->getTypeForMSH() - 1].numElem;
-        // Unique list of vertices
-        const int nVert = (*itElem)->getNumVertices();
-        for(int iVert = 0; iVert != nVert; ++iVert)
-          vertMap[(*itElem)->getVertex(iVert)] = 0;  // Unlabelled
-        // Maintain list of (base element) faces with only one bounding element.
-        for(int iFace = 0; iFace != ElemFaceTr<DIM, Elem>::numFaceT;
-            ++iFace) {
-          FaceT face = FaceTr<FaceT>::template getFace<Elem>(*itElem, iFace);
-          std::pair<typename BoFaceMap::iterator, bool> insBoFaceMap =
-            boFaceMap.insert(std::pair<FaceT, FaceData>
-                             (face, FaceData(*itElem, iFace, eVecIndex)));
-          if(!insBoFaceMap.second) {
-            // The face already exists and is therefore bounded on both sides
-            // by elements.  Not a boundary face so delete.
-            boFaceMap.erase(insBoFaceMap.first);
+    unsigned numElem[4];
+    numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+    entity->getNumMeshElements(numElem);
+    // Loop over all types of elements
+    int nType = entity->getNumElementTypes();
+    for(int iType = 0; iType != nType; ++iType) {
+      // Loop over all elements in a type
+      MElement *const *element = entity->getStartElementType(iType);
+      const int nElem = numElem[iType];
+      for(int iElem = 0; iElem != nElem; ++iElem) {
+        if(partition < 0 || element[iElem]->getPartition() == partition) {
+          // Unique list of elements
+          const unsigned eVecIndex = elemVec.size();
+          elemVec.push_back(ElemData(element[iElem]));
+          ++zoneElemConn[(element[iElem])->getTypeForMSH() - 1].numElem;
+          // Unique list of vertices
+          const int nVert = element[iElem]->getNumVertices();
+          for(int iVert = 0; iVert != nVert; ++iVert)
+            vertMap[element[iElem]->getVertex(iVert)] = 0;  // Unlabelled
+          // Maintain list of (base element) faces with only one bounding
+          // element.
+          const int nFace = DimTr<DIM>::getNumFace(element[iElem]);
+          for(int iFace = 0; iFace != nFace; ++iFace) {
+            FaceT face = DimTr<DIM>::getFace(element[iElem], iFace);
+            std::pair<typename BoFaceMap::iterator, bool> insBoFaceMap =
+              boFaceMap.insert(std::pair<FaceT, FaceData>
+                               (face, FaceData(element[iElem], iFace,
+                                               eVecIndex)));
+            if(!insBoFaceMap.second) {
+              // The face already exists and is therefore bounded on both sides
+              // by elements.  Not a boundary face so delete.
+              boFaceMap.erase(insBoFaceMap.first);
+            }
           }
         }
       }
diff --git a/Geo/MZone.h b/Geo/MZone.h
index 2d5aeebd50..7ea7bc1940 100644
--- a/Geo/MZone.h
+++ b/Geo/MZone.h
@@ -23,7 +23,6 @@
 #include "MEdgeHash.h"
 #include "MFaceHash.h"
 #include "CustomContainer.h"
-#include "ElementTraits.h"
 
 // #define HAVE_HASH_MAP
 
diff --git a/Geo/Makefile b/Geo/Makefile
index 9d272c1b8a..68c9c1103f 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -221,7 +221,7 @@ GModelIO_CGNS.o: GModelIO_CGNS.cpp GModel.h GVertex.h GEntity.h Range.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h MZone.h \
   MElement.h ../Common/GmshDefines.h MVertex.h MEdge.h MFace.h \
   MEdgeHash.h ../Common/Hash.h MFaceHash.h CustomContainer.h \
-  ElementTraits.h MZoneBoundary.h CGNSOptions.h
+  MZoneBoundary.h CGNSOptions.h
 GModelIO_MED.o: GModelIO_MED.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h SVector3.h \
   GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/Message.h
@@ -294,9 +294,9 @@ MZone.o: MZone.cpp MZone.h MElement.h ../Common/GmshDefines.h MVertex.h \
   SPoint3.h MEdge.h SVector3.h MFace.h GFace.h GEntity.h Range.h \
   SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h SPoint2.h \
   Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h MFaceHash.h \
-  CustomContainer.h ElementTraits.h
+  CustomContainer.h
 MZoneBoundary.o: MZoneBoundary.cpp MZoneBoundary.h MZone.h MElement.h \
   ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
   GFace.h GEntity.h Range.h SBoundingBox3d.h GPoint.h GEdgeLoop.h GEdge.h \
   GVertex.h SPoint2.h Pair.h GRegion.h MEdgeHash.h ../Common/Hash.h \
-  MFaceHash.h CustomContainer.h ElementTraits.h ../Common/Message.h
+  MFaceHash.h CustomContainer.h ../Common/Message.h
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 0566d60546..035da47004 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -391,4 +391,16 @@ HighOrder.o: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
   ../Numeric/NumericEmbedded.h ../Common/Context.h ../Geo/CGNSOptions.h \
   ../Mesh/PartitionOptions.h ../Common/GmshMatrix.h \
   ../Numeric/FunctionSpace.h
-Partition.o: Partition.cpp
+Partition.o: Partition.cpp ../Geo/MElement.h \
+  ../Common/GmshDefines.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/GEdge.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h \
+  ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
+  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
+  ../Geo/GRegion.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h Partition.h \
+  ../Common/Message.h PartitionOptions.h
diff --git a/Mesh/Partition.cpp b/Mesh/Partition.cpp
index 629d806ca4..03873774d4 100644
--- a/Mesh/Partition.cpp
+++ b/Mesh/Partition.cpp
@@ -7,7 +7,6 @@
 
 #if defined(HAVE_CHACO) || defined(HAVE_METIS)
 
-#include "ElementTraits.h"
 #include "GModel.h"
 #include "Partition.h"
 #include "PartitionObjects.h"
@@ -56,9 +55,9 @@ template <> struct LFaceTr<MFace>
  * Forward declarations
  *============================================================================*/
 
-template <unsigned DIM, unsigned N = DimTr<DIM>::numElemTypes>
+template <unsigned DIM>
 struct MakeGraphFromEntity;
-template <unsigned DIM, unsigned N = DimTr<DIM-1>::numElemTypes>
+template <unsigned DIM>
 struct MatchBoElemToGrVertex;
 template <unsigned DIM, typename EntIter, typename EntIterBE>
 void MakeGraphDIM(const EntIter begin, const EntIter end,
@@ -100,11 +99,21 @@ int PartitionMesh(GModel *const model, PartitionOptions &options)
     return 1;
   }
 
-  // Assign partitions to internal elements
+  // Count partition sizes and assign partitions to internal elements
+  std::vector<int> ssize(options.num_partitions, 0);
   const int n = graph.getNumVertex();
   for(int i = 0; i != n; ++i) {
+    ++ssize[graph.partition[i] - 1];
     graph.element[i]->setPartition(graph.partition[i]);
   }
+  int sMin = graph.getNumVertex();
+  int sMax = 0;
+  for(int i = 0; i != options.num_partitions; ++i) {
+    sMin = std::min(sMin, ssize[i]);
+    sMax = std::max(sMax, ssize[i]);
+  }
+  model->setMinPartitionSize(sMin);
+  model->setMaxPartitionSize(sMax);
 
   // Assign partitions to boundary elements
   const int nb = boElemGrVec.size();
@@ -148,6 +157,22 @@ int PartitionGraph(Graph &graph, PartitionOptions &options)
           options.vmax = 2*(1 << options.ndims);
 	}
       }
+      // Ensure num_partitions reflects values used by Chaco
+      switch(options.architecture) {
+      case 0:
+        options.num_partitions = 1 << options.ndims_tot;
+        break;
+      case 1:
+        options.num_partitions = options.mesh_dims[0];
+        break;
+      case 2:
+        options.num_partitions = options.mesh_dims[0]*options.mesh_dims[1];
+        break;
+      case 3:
+        options.num_partitions =
+          options.mesh_dims[0]*options.mesh_dims[1]*options.mesh_dims[2];
+        break;
+      }
       try {
         const int iSec = 0;
         ier = interface
@@ -168,8 +193,8 @@ int PartitionGraph(Graph &graph, PartitionOptions &options)
         ier = 2;
       }
       if(!ier) graph.short2int();
-      return ier;
     }
+    break;
 #endif
   case 2:  // Metis
 #ifdef HAVE_METIS
@@ -219,15 +244,14 @@ int PartitionGraph(Graph &graph, PartitionOptions &options)
         ier = 2;
       }
       // Number partitions from 1
-      {
+      if(!ier) {
         int *p = &graph.partition[0];  //**Sections
         for(int n = graph.getNumVertex(); n--;) ++(*p++);
       }
-      return ier;
     }
+    break;
 #endif
   }
-
   return ier;
 
 }
@@ -279,28 +303,11 @@ int MakeGraph(GModel *const model, Graph &graph, BoElemGrVec *const boElemGrVec)
 
 //--Get the dimension of the mesh and count the numbers of elements
 
-      int numElem[4];
-      int meshDim = 3;
-      for(int i = 0; i != 4; ++i) numElem[i] = 0;
-      for(GModel::riter it = model->firstRegion(); it != model->lastRegion();
-          ++it) {
-        numElem[ElemTypeTetra] += (*it)->tetrahedra.size();
-        numElem[ElemTypeHexa] += (*it)->hexahedra.size();
-        numElem[ElemTypePrism] += (*it)->prisms.size();
-        numElem[ElemTypePyramid] += (*it)->pyramids.size();
-      }
-      if(numElem[ElemTypeTetra] + numElem[ElemTypeHexa] +
-         numElem[ElemTypePrism] + numElem[ElemTypePyramid] == 0) {
-        for(GModel::fiter it = model->firstFace(); it != model->lastFace();
-            ++it) {
-          numElem[ElemTypeTri] += (*it)->triangles.size();
-          numElem[ElemTypeQuad] += (*it)->quadrangles.size();
-        }
-        if(numElem[ElemTypeTri] + numElem[ElemTypeQuad] == 0) {
-          Msg::Error("No mesh elements were found");
-          return 1;
-        }
-        meshDim = 2;
+      unsigned numElem[4];
+      const int meshDim = model->getNumMeshElements(numElem);
+      if(meshDim < 2) {
+        Msg::Error("No mesh elements were found");
+        return 1;
       }
 
 //--Make the graph
@@ -479,8 +486,6 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
 {
 
   typedef typename DimTr<DIM>::FaceT FaceT;
-  typedef typename DimTr<DIM>::EntityT Ent;
-  typedef typename DimTr<DIM-1>::EntityT EntBE;
   typedef typename LFaceTr<FaceT>::FaceMap FaceMap;
 
   graph.markSection();
@@ -489,8 +494,7 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
   GrVertexMap grVertMap;
 
   for(EntIter entIt = begin; entIt != end; ++entIt) {
-    MakeGraphFromEntity<DIM>::eval
-       (static_cast<Ent*>(*entIt), faceMap, grVertMap, graph);
+    MakeGraphFromEntity<DIM>::eval(*entIt, faceMap, grVertMap, graph);
   }
 
   // Write any graph vertices remaining in 'grVertMap'.  These are boundary
@@ -507,7 +511,7 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
     boElemGrVec->reserve(faceMap.size());
     for(EntIterBE entIt = beginBE; entIt != endBE; ++entIt) {
       MatchBoElemToGrVertex<DIM>::eval
-        (static_cast<EntBE*>(*entIt), faceMap, grVertMap, graph, *boElemGrVec);
+        (*entIt, faceMap, grVertMap, graph, *boElemGrVec);
     }
   }
 
@@ -523,132 +527,66 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
  *
  *   Helps generate a graph - operates on a single entity
  *
- * Notes
- * =====
- *
- *   - template meta-programming is used to iterate over the various types of
- *     elements in an entity
- *
  ******************************************************************************/
 
-//--Entry point
-
-template <unsigned DIM, unsigned N>
+template <unsigned DIM>
 struct MakeGraphFromEntity
 {
-
-  typedef typename DimTr<DIM>::EntityT Ent;  // The type of entity
   typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
   typedef typename LFaceTr<FaceT>::FaceMap FaceMap;  // The corresponding map
-  typedef typename DimElemTr<DIM, N>::Elem Elem; // The type of primary element
-  typedef typename DimElemTr<DIM, N>::ConstElementIterator ConstElementIterator;
-      
-  static void eval(const Ent *const entity, FaceMap &faceMap,
-                   GrVertexMap &grVertMap, Graph &graph)
-  {
-    // Loop over all elements in the entity
-    ConstElementIterator elemEnd = DimElemTr<DIM, N>::end(entity);
-    for(ConstElementIterator elemIt = DimElemTr<DIM, N>::begin(entity);
-        elemIt != elemEnd; ++elemIt) {
-      // Insert this element into the map of graph vertices
-      std::pair<typename GrVertexMap::iterator, bool> insGrVertMap =
-        grVertMap.insert
-        (std::pair<MElement*, GrVertex>
-         (*elemIt, GrVertex(graph.getNextIndex(),
-                            ElemFaceTr<DIM, Elem>::numFaceT)));
-      // Loop over all faces in the element
-      for(int iFace = 0; iFace != ElemFaceTr<DIM, Elem>::numFaceT; ++iFace) {
-        FaceT face = FaceTr<FaceT>::template getFace<Elem>(*elemIt, iFace);
-        std::pair<typename FaceMap::iterator, bool> insFaceMap =
-          faceMap.insert(std::pair<FaceT, MElement*>(face, *elemIt));
-        if(!insFaceMap.second) {
-          // The face already exists
-          typename GrVertexMap::iterator grVertIt2 =
-            grVertMap.find(insFaceMap.first->second);
-                                        // Iterator to the second graph vertex
-                                        // that connects to this face
-          // Update edges in both graph vertices
-          grVertIt2->second.add(insGrVertMap.first->second.index);
-          insGrVertMap.first->second.add(grVertIt2->second.index);
-          if(grVertIt2->second.complete()) {
-            // This graph vertex has complete connectivity.  Write and delete.
-            graph.add(grVertIt2);
-            grVertMap.erase(grVertIt2);
-          }
-          // The face is no longer required
-          faceMap.erase(insFaceMap.first);
-        }
-      }
-      if(insGrVertMap.first->second.complete()) {
-        // This graph vertex already has complete connectivity.  Write and
-        // delete.
-        graph.add(insGrVertMap.first);
-        grVertMap.erase(insGrVertMap.first);
-      }
-    }
-    // Next element type in the entity
-    MakeGraphFromEntity<DIM, N-1>::eval(entity, faceMap, grVertMap, graph);
-  }
 
-};
-
-//--Terminate loop when N = 1
-
-template <unsigned DIM>
-struct MakeGraphFromEntity<DIM, 1>
-{
-
-  typedef typename DimTr<DIM>::EntityT Ent;  // The type of entity
-  typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
-  typedef typename LFaceTr<FaceT>::FaceMap FaceMap;  // The corresponding map
-  typedef typename DimElemTr<DIM, 1>::Elem Elem; // The type of primary element
-  typedef typename DimElemTr<DIM, 1>::ConstElementIterator ConstElementIterator;
-      
-  static void eval(const Ent *const entity, FaceMap &faceMap,
+  static void eval(const GEntity *const entity, FaceMap &faceMap,
                    GrVertexMap &grVertMap, Graph &graph)
   {
-    // Loop over all elements in the entity
-    ConstElementIterator elemEnd = DimElemTr<DIM, 1>::end(entity);
-    for(ConstElementIterator elemIt = DimElemTr<DIM, 1>::begin(entity);
-        elemIt != elemEnd; ++elemIt) {
-      // Insert this element into the map of graph vertices
-      std::pair<typename GrVertexMap::iterator, bool> insGrVertMap =
-        grVertMap.insert
-        (std::pair<MElement*, GrVertex>
-         (*elemIt, GrVertex(graph.getNextIndex(),
-                            ElemFaceTr<DIM, Elem>::numFaceT)));
-      // Loop over all faces in the element
-      for(int iFace = 0; iFace != ElemFaceTr<DIM, Elem>::numFaceT; ++iFace) {
-        FaceT face = FaceTr<FaceT>::template getFace<Elem>(*elemIt, iFace);
-        std::pair<typename FaceMap::iterator, bool> insFaceMap =
-          faceMap.insert(std::pair<FaceT, MElement*>(face, *elemIt));
-        if(!insFaceMap.second) {
-          // The face already exists
-          typename GrVertexMap::iterator grVertIt2 =
-            grVertMap.find(insFaceMap.first->second);
+    unsigned numElem[4];
+    numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+    entity->getNumMeshElements(numElem);
+    // Loop over all types of elements
+    int nType = entity->getNumElementTypes();
+    for(int iType = 0; iType != nType; ++iType) {
+      // Loop over all elements in a type
+      MElement *const *element = entity->getStartElementType(iType);
+      const int nElem = numElem[iType];
+      for(int iElem = 0; iElem != nElem; ++iElem) {
+        const int nFace = DimTr<DIM>::getNumFace(element[iElem]);
+        // Insert this element into the map of graph vertices
+        std::pair<typename GrVertexMap::iterator, bool> insGrVertMap =
+          grVertMap.insert
+          (std::pair<MElement*, GrVertex>
+           (element[iElem], GrVertex(graph.getNextIndex(), nFace)));
+        // Loop over all faces in the element
+        for(int iFace = 0; iFace != nFace; ++iFace) {
+          FaceT face = DimTr<DIM>::getFace(element[iElem], iFace);
+          std::pair<typename FaceMap::iterator, bool> insFaceMap =
+             faceMap.insert(std::pair<FaceT, MElement*>(face, element[iElem]));
+          if(!insFaceMap.second) {
+             // The face already exists
+             typename GrVertexMap::iterator grVertIt2 =
+                grVertMap.find(insFaceMap.first->second);
                                         // Iterator to the second graph vertex
                                         // that connects to this face
-          // Update edges in both graph vertices
-          grVertIt2->second.add(insGrVertMap.first->second.index);
-          insGrVertMap.first->second.add(grVertIt2->second.index);
-          if(grVertIt2->second.complete()) {
-            // This graph vertex has complete connectivity.  Write and delete.
-            graph.add(grVertIt2);
-            grVertMap.erase(grVertIt2);
+             // Update edges in both graph vertices
+             grVertIt2->second.add(insGrVertMap.first->second.index);
+             insGrVertMap.first->second.add(grVertIt2->second.index);
+             if(grVertIt2->second.complete()) {
+                // This graph vertex has complete connectivity.  Write and
+                // delete.
+                graph.add(grVertIt2);
+                grVertMap.erase(grVertIt2);
+             }
+             // The face is no longer required
+             faceMap.erase(insFaceMap.first);
           }
-          // The face is no longer required
-          faceMap.erase(insFaceMap.first);
         }
-      }
-      if(insGrVertMap.first->second.complete()) {
-        // This graph vertex already has complete connectivity.  Write and
-        // delete.
-        graph.add(insGrVertMap.first);
-        grVertMap.erase(insGrVertMap.first);
+        if(insGrVertMap.first->second.complete()) {
+           // This graph vertex already has complete connectivity.  Write and
+           // delete.
+           graph.add(insGrVertMap.first);
+           grVertMap.erase(insGrVertMap.first);
+        }
       }
     }
   }
-
 };
 
 
@@ -664,75 +602,38 @@ struct MakeGraphFromEntity<DIM, 1>
  *
  ******************************************************************************/
 
-//--Entry point
-
-template <unsigned DIM, unsigned N>
-struct MatchBoElemToGrVertex
-{
-
-  typedef typename DimTr<DIM-1>::EntityT Ent;  // The type of entity
-  typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
-  typedef typename LFaceTr<FaceT>::FaceMap FaceMap;  // The corresponding map
-  typedef typename DimElemTr<DIM-1, N>::Elem Elem; // The type of primary Elem.
-  typedef typename DimElemTr<DIM-1, N>::ConstElementIterator
-    ConstElementIterator;
-      
-  static void eval(const Ent *const entity, const FaceMap &faceMap,
-                   const GrVertexMap &grVertMap, const Graph &graph,
-                   std::vector<BoElemGr> &boElemGrVec)
-  {
-    // Loop over all elements in the entity
-    const ConstElementIterator elemEnd = DimElemTr<DIM-1, N>::end(entity);
-    for(ConstElementIterator elemIt = DimElemTr<DIM-1, N>::begin(entity);
-        elemIt != elemEnd; ++elemIt) {
-      FaceT face = FaceTr<FaceT>::template getFace<Elem>(*elemIt, 0);
-      const typename FaceMap::const_iterator faceMapIt = faceMap.find(face);
-      if(faceMapIt != faceMap.end()) {
-        const typename GrVertexMap::const_iterator grVertMapIt =
-          grVertMap.find(faceMapIt->second);
-        boElemGrVec.push_back
-          (BoElemGr(*elemIt, graph.convertC2W(grVertMapIt->second.index) - 1));
-      }
-    }
-    // Next element type in the entity
-    MatchBoElemToGrVertex<DIM, N-1>::eval(entity, faceMap, grVertMap, graph,
-                                          boElemGrVec);
-  }
-
-};
-
-//--Terminate loop when N = 1
-
 template <unsigned DIM>
-struct MatchBoElemToGrVertex<DIM, 1>
+struct MatchBoElemToGrVertex
 {
-
-  typedef typename DimTr<DIM-1>::EntityT Ent;  // The type of entity
   typedef typename DimTr<DIM>::FaceT FaceT;  // The type/dimension of face
   typedef typename LFaceTr<FaceT>::FaceMap FaceMap;  // The corresponding map
-  typedef typename DimElemTr<DIM-1, 1>::Elem Elem; // The type of primary Elem.
-  typedef typename DimElemTr<DIM-1, 1>::ConstElementIterator
-    ConstElementIterator;
       
-  static void eval(const Ent *const entity, const FaceMap &faceMap,
+  static void eval(const GEntity *const entity, const FaceMap &faceMap,
                    const GrVertexMap &grVertMap, const Graph &graph,
                    std::vector<BoElemGr> &boElemGrVec)
   {
-    // Loop over all elements in the entity
-    const ConstElementIterator elemEnd = DimElemTr<DIM-1, 1>::end(entity);
-    for(ConstElementIterator elemIt = DimElemTr<DIM-1, 1>::begin(entity);
-        elemIt != elemEnd; ++elemIt) {
-      FaceT face = FaceTr<FaceT>::template getFace<Elem>(*elemIt, 0);
-      const typename FaceMap::const_iterator faceMapIt = faceMap.find(face);
-      if(faceMapIt != faceMap.end()) {
-        const typename GrVertexMap::const_iterator grVertMapIt =
-          grVertMap.find(faceMapIt->second);
-        boElemGrVec.push_back
-          (BoElemGr(*elemIt, graph.convertC2W(grVertMapIt->second.index) - 1));
+    unsigned numElem[4];
+    numElem[0] = 0; numElem[1] = 0; numElem[2] = 0; numElem[3] = 0;
+    entity->getNumMeshElements(numElem);
+    // Loop over all types of elements
+    int nType = entity->getNumElementTypes();
+    for(int iType = 0; iType != nType; ++iType) {
+      // Loop over all elements in a type
+      MElement *const *element = entity->getStartElementType(iType);
+      const int nElem = numElem[iType];
+      for(int iElem = 0; iElem != nElem; ++iElem) {
+        FaceT face = DimTr<DIM>::getFace(element[iElem], 0);
+        const typename FaceMap::const_iterator faceMapIt = faceMap.find(face);
+        if(faceMapIt != faceMap.end()) {
+          const typename GrVertexMap::const_iterator grVertMapIt =
+            grVertMap.find(faceMapIt->second);
+          boElemGrVec.push_back
+            (BoElemGr(element[iElem],
+                      graph.convertC2W(grVertMapIt->second.index) - 1));
+        }
       }
     }
   }
-
 };
 
 
diff --git a/contrib/Chaco/main/Gmsh_printf.cpp b/contrib/Chaco/main/Gmsh_printf.cpp
index 167775c5e4..516ae649e6 100644
--- a/contrib/Chaco/main/Gmsh_printf.cpp
+++ b/contrib/Chaco/main/Gmsh_printf.cpp
@@ -3,6 +3,7 @@
 // See the LICENSE.txt file for license information. Please report all
 // bugs and problems to <gmsh@geuz.org>.
 
+#include <cstring>
 #include "Message.h"
 
 // Overload the printf statements in Chaco to write using Msg::Direct in gmsh
-- 
GitLab