Skip to content
Snippets Groups Projects
Select Git revision
  • 6371397b414e875a0bd2f77527d35061b7833d08
  • master default
  • library-names
  • fix_script_header
  • fix_libdir
  • fix_cmake_hdf5
  • partition
  • cgnsUnstructured
  • partitioning
  • HighOrderBLCurving
  • gmsh_3_0_5
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
30 results

Callbacks.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    meshPartition.cpp 90.62 KiB
    // Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <gmsh@onelab.info>.
    //
    // Contributed by Anthony Royer.
    
    #include <vector>
    #include <set>
    #include <sstream>
    #include <algorithm>
    #include <ctime>
    #include <limits>
    #include <stack>
    #include <cstdlib>
    #include <map>
    #include <utility>
    #include "GmshConfig.h"
    #include "GmshMessage.h"
    #include "GModel.h"
    
    // TODO C++11 remove this nasty stuff
    #if __cplusplus >= 201103L
    #include <unordered_map>
    #define hashmap std::unordered_map
    #define hashmapface std::unordered_map\
      <MFace, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
       Hash_Face, Equal_Face>
    #define hashmapedge std::unordered_map\
      <MEdge, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
       Hash_Edge, Equal_Edge>
    #define hashmapvertex std::unordered_map\
      <MVertex*, std::vector<std::pair<MElement*, std::vector<unsigned int> > > >
    #else
    #define hashmap std::map
    #define hashmapface std::map\
      <MFace, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
       Less_Face>
    #define hashmapedge std::map\
      <MEdge, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
       Less_Edge>
    #define hashmapvertex std::map\
      <MVertex*, std::vector<std::pair<MElement*, std::vector<unsigned int> > > >
    #endif
    
    #if defined(HAVE_METIS)
    
    #include "OS.h"
    #include "Context.h"
    #include "partitionRegion.h"
    #include "partitionFace.h"
    #include "partitionEdge.h"
    #include "partitionVertex.h"
    #include "ghostRegion.h"
    #include "ghostFace.h"
    #include "ghostEdge.h"
    #include "MFaceHash.h"
    #include "MEdgeHash.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MTetrahedron.h"
    #include "MHexahedron.h"
    #include "MPrism.h"
    #include "MPyramid.h"
    #include "MTrihedron.h"
    #include "MElementCut.h"
    #include "MPoint.h"
    
    extern "C" {
    #include <metis.h>
    }
    
    // Graph of the mesh for partitioning purposes.
    class Graph
    {
     private:
      // The GModel
      GModel *_model;
      // The number of partitions
      unsigned int _nparts;
      // The number of elements
      size_t _ne;
      // The number of nodes
      size_t _nn;
      // The dimension of the mesh
      unsigned int _dim;
      // The list of nodes belonging to the ith element of the mesh are stored in
      // consecutive locations of eind starting at position eptr[i] up to (but not
      // including) position eptr[i+1]. The size of the eind array is of size equal
      // to the sum of the number of nodes in all the elements of the mesh.
      std::vector<size_t> _eind;
      // The size of the eptr array is n + 1, where n is the number of elements in
      // the mesh.
      std::vector<size_t> _eptr;
      // The metis graph structure
      std::vector<size_t> _xadj, _adjncy;
      // Elements corresponding to each graph elements in eptr
      std::vector<MElement*> _element;
      // Vertices corresponding to each graph vertices in eptr
      std::vector<int> _vertex;
      // The width associated to each elements of eptr
      std::vector<unsigned int> _vwgt;
      // The partitions output from the partitioner
      std::vector<unsigned int> _partition;
     public:
      Graph(GModel * const model)
        : _model(model), _nparts(0), _ne(0), _nn(0), _dim(0), _eind(), _eptr(),
        _xadj(), _adjncy(), _element(), _vertex(), _vwgt(), _partition()
      {
      }
      void fillDefaultWeights()
      {
        if(CTX::instance()->mesh.partitionTriWeight == 1 &&
           CTX::instance()->mesh.partitionQuaWeight == 1 &&
           CTX::instance()->mesh.partitionTetWeight == 1 &&
           CTX::instance()->mesh.partitionPyrWeight == 1 &&
           CTX::instance()->mesh.partitionPriWeight == 1 &&
           CTX::instance()->mesh.partitionHexWeight == 1) return;
    
        _vwgt.resize(_ne);
        for(unsigned int i = 0; i < _ne; i++){
          if(!_element[i]){
            _vwgt[i] = 1;
            continue;
          }
          switch (_element[i]->getType()){
          case TYPE_TRI: _vwgt[i] = CTX::instance()->mesh.partitionTriWeight; break;
          case TYPE_QUA: _vwgt[i] = CTX::instance()->mesh.partitionQuaWeight; break;
          case TYPE_TET: _vwgt[i] = CTX::instance()->mesh.partitionTetWeight; break;
          case TYPE_PYR: _vwgt[i] = CTX::instance()->mesh.partitionPyrWeight; break;
          case TYPE_PRI: _vwgt[i] = CTX::instance()->mesh.partitionPriWeight; break;
          case TYPE_HEX: _vwgt[i] = CTX::instance()->mesh.partitionHexWeight; break;
          default: _vwgt[i] = 1; break;
          }
        }
      }
      ~Graph()
      {
        clear();
      }
      unsigned int nparts() const { return _nparts; };
      size_t ne() const { return _ne; };
      size_t nn() const { return _nn; };
      unsigned int dim() const { return _dim; };
      size_t eind(size_t i) const { return _eind[i]; };
      size_t eptr(size_t i) const { return _eptr[i]; };
      size_t xadj(size_t i) const { return _xadj[i]; };
      std::vector<size_t> &xadj() { return _xadj; };
      size_t adjncy(size_t i) const { return _adjncy[i]; };
      std::vector<size_t> &adjncy() { return _adjncy; };
      MElement* element(size_t i) const { return _element[i]; };
      int vertex(size_t i) const { return _vertex[i]; };
      std::vector<unsigned int> &vwgt() { return _vwgt; };
      unsigned int partition(unsigned int i) const { return _partition[i]; };
      std::vector<unsigned int> &partition() { return _partition; };
      size_t numNodes() const { return _ne; };
      size_t numEdges() const { return _xadj[_ne]/2; };
      void nparts(unsigned int nparts) { _nparts = nparts; };
      void ne(size_t ne) { _ne = ne; };
      void nn(size_t nn) { _nn = nn; };
      void dim(unsigned int dim) { _dim = dim; };
      void eindResize(size_t size) { _eind.resize(size,0); }
      void eind(size_t i, size_t eind) { _eind[i] = eind; };
      void eptrResize(size_t size) { _eptr.resize(size, 0); }
      void eptr(size_t i, size_t eptr) { _eptr[i] = eptr; };
      void elementResize(size_t size){ _element.resize(size, 0); }
      void element(size_t i, MElement* element) { _element[i] = element; };
      void vertexResize(size_t size) { _vertex.resize(size, -1); }
      void adjncy(size_t i, size_t adjncy) { _adjncy[i] = adjncy; };
      void vertex(size_t i, int vertex) { _vertex[i] = vertex; };
      void vwgt(std::vector<unsigned int> &vwgt) { std::swap(_vwgt, vwgt); };
      void partition(std::vector<unsigned int> &partition) { std::swap(_partition, partition); };
      void clear()
      {
        _eind.clear();
        _eptr.clear();
        _xadj.clear();
        _adjncy.clear();
        _element.clear();
        _vertex.clear();
        _vwgt.clear();
        _partition.clear();
      }
      void clearDualGraph()
      {
        _xadj.clear();
        _adjncy.clear();
      }
      void eraseVertex()
      {
        for(size_t i = 0; i < _vertex.size(); ++i) _vertex[i] = -1;
      }
    
      std::vector< std::set<MElement*> > getBoundaryElements(size_t size = 0)
      {
        std::vector< std::set<MElement*> > elements
          ((size ? size : _nparts), std::set<MElement*>());
        for(size_t i = 0; i < _ne; ++i){
          for(size_t j = _xadj[i]; j < _xadj[i+1]; ++j){
            if(_partition[i] != _partition[_adjncy[j]]){
              if(_element[i]->getDim() == (int)_dim){
                elements[_partition[i]].insert(_element[i]);
              }
            }
          }
        }
    
        return elements;
      }
    
      void assignGhostCells()
      {
        std::vector<GEntity*> ghostEntities(_nparts, (GEntity*)NULL);
        int elementaryNumber = _model->getMaxElementaryNumber(_dim);
    
        for(unsigned int i = 1; i <= _nparts; ++i){
          switch (_dim) {
            case 1:
              ghostEntities[i - 1] = new ghostEdge(_model, ++elementaryNumber, i);
              _model->add(static_cast<ghostEdge*>(ghostEntities[i - 1]));
              break;
            case 2:
              ghostEntities[i - 1] = new ghostFace(_model, ++elementaryNumber, i);
              _model->add(static_cast<ghostFace*>(ghostEntities[i - 1]));
              break;
            case 3:
              ghostEntities[i - 1] = new ghostRegion(_model, ++elementaryNumber, i);
              _model->add(static_cast<ghostRegion*>(ghostEntities[i - 1]));
              break;
            default:
              break;
          }
        }
    
        for(size_t i = 0; i < _ne; ++i){
          std::set<short> ghostCellsPartition;
          for(size_t j = _xadj[i]; j < _xadj[i+1]; j++){
            if(_partition[i] != _partition[_adjncy[j]] &&
               ghostCellsPartition.find(_partition[_adjncy[j]]) == ghostCellsPartition.end()){
              if(_element[i]->getDim() == (int)_dim){
                switch (_dim) {
                  case 1:
                    static_cast<ghostEdge*>(ghostEntities[_partition[_adjncy[j]]])->
                    addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
                    break;
                  case 2:
                    static_cast<ghostFace*>(ghostEntities[_partition[_adjncy[j]]])->
                    addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
                    break;
                  case 3:
                    static_cast<ghostRegion*>(ghostEntities[_partition[_adjncy[j]]])->
                    addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
                    break;
                  default:
                    break;
                }
                ghostCellsPartition.insert(_partition[_adjncy[j]]);
              }
            }
          }
        }
      }
    
      void createDualGraph(bool connectedAll)
      {
        std::vector<size_t> nptr(_nn+1, 0);
        std::vector<size_t> nind(_eptr[_ne], 0);
    
        for(size_t i = 0; i < _ne; ++i){
          for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
            nptr[_eind[j]]++;
          }
        }
    
        for(size_t i = 1; i < _nn; ++i) nptr[i] += nptr[i-1];
        for(size_t i = _nn; i > 0; --i) nptr[i] = nptr[i-1];
        nptr[0] = 0;
    
        for(size_t i = 0; i < _ne; ++i){
          for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
            nind[nptr[_eind[j]]++] = i;
          }
        }
    
        for(unsigned int i = _nn; i > 0; --i) nptr[i] = nptr[i-1];
        nptr[0] = 0;
    
        _xadj.resize(_ne+1, 0);
        std::vector<size_t> nbrs(_ne, 0);
        std::vector<unsigned int> marker(_ne, 0);
    
        for(size_t i = 0; i < _ne; ++i){
          size_t l = 0;
          for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
            for(size_t k = nptr[_eind[j]]; k < nptr[_eind[j]+1]; ++k){
              if(nind[k] != i){
                if(marker[nind[k]] == 0) nbrs[l++] = nind[k];
                marker[nind[k]]++;
              }
            }
          }
    
          unsigned int nbrsNeighbors = 0;
          for(size_t j = 0; j < l; j++){
            if(marker[nbrs[j]] >=
               (connectedAll ? 1 :
                _element[i]->numCommonNodesInDualGraph(_element[nbrs[j]])))
              nbrsNeighbors++;
            marker[nbrs[j]] = 0;
            nbrs[j] = 0;
          }
    
          _xadj[i] = nbrsNeighbors;
        }
    
        for(size_t i = 1; i < _ne; ++i) _xadj[i] = _xadj[i] + _xadj[i-1];
        for(size_t i = _ne; i > 0; --i) _xadj[i] = _xadj[i-1];
        _xadj[0] = 0;
    
        _adjncy.resize(_xadj[_ne], 0);
    
        for(size_t i = 0; i < _ne; ++i){
          size_t l = 0;
          for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
            for(size_t k = nptr[_eind[j]]; k < nptr[_eind[j]+1]; ++k){
              if(nind[k] != i){
                if (marker[nind[k]] == 0) nbrs[l++] = nind[k];
                marker[nind[k]]++;
              }
            }
          }
    
          for(size_t j = 0; j < l; ++j){
            if(marker[nbrs[j]] >=
               (connectedAll ? 1 :
                _element[i]->numCommonNodesInDualGraph(_element[nbrs[j]]))){
                 _adjncy[_xadj[i]] = nbrs[j];
                 _xadj[i] = _xadj[i]+1;
            }
            marker[nbrs[j]] = 0;
            nbrs[j] = 0;
          }
        }
        for(size_t i = _ne; i > 0; --i) _xadj[i] = _xadj[i-1];
        _xadj[0] = 0;
      }
    };
    
    template <class ITERATOR>
    static void fillElementsToNodesMap(Graph &graph, const GEntity *const entity,
                                       size_t &eptrIndex, size_t &eindIndex, size_t &numVertex,
                                       ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        const int numVertices = (*it)->getNumPrimaryVertices();
        graph.element(eptrIndex++, *it);
        graph.eptr(eptrIndex, graph.eptr(eptrIndex-1) + numVertices);
        for(int i = 0; i < numVertices; i++){
          if(graph.vertex((*it)->getVertex(i)->getNum()-1) == -1){
            graph.vertex((*it)->getVertex(i)->getNum()-1, numVertex++);
          }
          graph.eind(eindIndex, graph.vertex((*it)->getVertex(i)->getNum()-1));
          eindIndex++;
        }
      }
    }
    
    static size_t getSizeOfEind(const GModel *const model)
    {
      size_t size = 0;
      // Loop over regions
      for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
        size += 4*(*it)->tetrahedra.size();
        size += 8*(*it)->hexahedra.size();
        size += 6*(*it)->prisms.size();
        size += 5*(*it)->pyramids.size();
        size += 4*(*it)->trihedra.size();
      }
    
      // Loop over faces
      for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
        size += 3*(*it)->triangles.size();
        size += 4*(*it)->quadrangles.size();
      }
    
      // Loop over edges
      for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
        size += 2*(*it)->lines.size();
      }
    
      // Loop over vertices
      for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
        size += 1*(*it)->points.size();
      }
    
      return size;
    }
    
    // Creates a mesh data structure used by Metis routines. Returns: 0 = success, 1
    // = no elements found, 2 = error.
    static int MakeGraph(GModel *const model, Graph &graph, int selectDim)
    {
      size_t eindSize = 0;
      if(selectDim < 0){
        graph.ne(model->getNumMeshElements());
        graph.nn(model->getNumMeshVertices());
        graph.dim(model->getDim());
        graph.elementResize(graph.ne());
        graph.vertexResize(model->getMaxVertexNumber());
        graph.eptrResize(graph.ne()+1);
        graph.eptr(0,0);
        eindSize = getSizeOfEind(model);
        graph.eindResize(eindSize);
      }
      else{
        GModel* tmp = new GModel();
        std::vector<GEntity*> entities;
        model->getEntities(entities);
    
        std::set<MVertex*> vertices;
        for(unsigned int  i = 0; i < entities.size(); i++){
          if(entities[i]->dim() == selectDim){
            switch (entities[i]->dim()) {
              case 3:
                tmp->add(static_cast<GRegion*>(entities[i]));
                break;
              case 2:
                tmp->add(static_cast<GFace*>(entities[i]));
                break;
              case 1:
                tmp->add(static_cast<GEdge*>(entities[i]));
                break;
              case 0:
                tmp->add(static_cast<GVertex*>(entities[i]));
                break;
              default:
                break;
            }
            for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
              for(unsigned int k = 0; k < entities[i]->getMeshElement(j)->getNumVertices() ; k++){
                vertices.insert(entities[i]->getMeshElement(j)->getVertex(k));
              }
            }
          }
        }
    
        graph.ne(tmp->getNumMeshElements());
        graph.nn(vertices.size());
        graph.dim(tmp->getDim());
        graph.elementResize(graph.ne());
        graph.vertexResize(model->getMaxVertexNumber());
        graph.eptrResize(graph.ne()+1);
        graph.eptr(0,0);
        eindSize = getSizeOfEind(tmp);
        graph.eindResize(eindSize);
    
        tmp->remove();
        delete tmp;
      }
    
      size_t eptrIndex = 0;
      size_t eindIndex = 0;
      size_t numVertex = 0;
    
      if(graph.ne() == 0){
        Msg::Error("No mesh elements were found");
        return 1;
      }
      if(graph.dim() == 0){
        Msg::Error("Cannot partition a point");
        return 1;
      }
    
      // Loop over regions
      if(selectDim < 0 || selectDim == 3){
        for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
          const GRegion *r = *it;
          fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
                                 r->tetrahedra.begin(), r->tetrahedra.end());
          fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
                                 r->hexahedra.begin(), r->hexahedra.end());
          fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
                                 r->prisms.begin(), r->prisms.end());
          fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
                                 r->pyramids.begin(), r->pyramids.end());
          fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
                                 r->trihedra.begin(), r->trihedra.end());
        }
      }
    
      // Loop over faces
      if(selectDim < 0 || selectDim == 2){
        for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
          const GFace *f = *it;
    
          fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
                                 f->triangles.begin(), f->triangles.end());
          fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
                                 f->quadrangles.begin(), f->quadrangles.end());
        }
      }
    
      // Loop over edges
      if(selectDim < 0 || selectDim == 1){
        for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
          const GEdge *e = *it;
          fillElementsToNodesMap(graph, e, eptrIndex, eindIndex, numVertex,
                                 e->lines.begin(), e->lines.end());
        }
      }
    
      // Loop over vertices
      if(selectDim < 0 || selectDim == 0){
        for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
          GVertex *v = *it;
    
          fillElementsToNodesMap(graph, v, eptrIndex, eindIndex, numVertex,
                                 v->points.begin(), v->points.end());
        }
      }
    
      return 0;
    }
    
    // Partition a graph created by MakeGraph using Metis library. Returns: 0 =
    // success, 1 = error, 2 = exception thrown.
    static int PartitionGraph(Graph &graph)
    {
    #ifdef HAVE_METIS
      Msg::Info("Running METIS graph partitioner");
    
      try {
        int metisOptions[METIS_NOPTIONS];
        METIS_SetDefaultOptions((idx_t *)metisOptions);
    
        switch(CTX::instance()->mesh.metisAlgorithm){
          case 1: // Recursive
            metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_RB;
            break;
          case 2: // K-way
            metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
            break;
          default:
            Msg::Info("Unknown partitioning algorithm");
            break;
        }
    
        switch(CTX::instance()->mesh.metisEdgeMatching){
          case 1: // Random matching
            metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_RM;
            break;
          case 2: // Sorted heavy-edge matching
            metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM;
            break;
          default:
            Msg::Info("Unknown partitioning edge matching");
            break;
        }
    
        switch(CTX::instance()->mesh.metisRefinementAlgorithm){
          case 1: // FM-based cut refinement
            metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_FM;
            break;
          case 2: // Greedy boundary refinement
            metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY;
            break;
          case 3: // Two-sided node FM refinement
            metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP2SIDED;
            break;
          case 4: // One-sided node FM refinement
            metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP1SIDED;
            break;
          default:
            Msg::Info("Unknown partitioning refinement algorithm");
            break;
        }
    
        // C numbering
        metisOptions[METIS_OPTION_NUMBERING] = 0;
        // Specifies the type of objective
        metisOptions[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
        // Forces contiguous partitions.
        //metisOptions[METIS_OPTION_CONTIG] = 1;
    
        int objval;
        std::vector<unsigned int> epart(graph.ne());
        const unsigned int ne = graph.ne();
        const int numPart = graph.nparts();
        int ncon = 1;
    
        graph.fillDefaultWeights();
    
        int metisError = 0;
        graph.createDualGraph(false);
    
        if (metisOptions[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY){
          metisError = METIS_PartGraphKway
            ((idx_t *)&ne, (idx_t *)&ncon, (idx_t *)graph.xadj().data(),
             (idx_t *)graph.adjncy().data(), (idx_t *)graph.vwgt().data(),
             (idx_t *)0, 0, (idx_t *)&numPart, 0, 0,
             (idx_t *)metisOptions, (idx_t *)&objval, (idx_t *)epart.data());
        }
        else{
          metisError = METIS_PartGraphRecursive
            ((idx_t *)&ne, (idx_t *)&ncon, (idx_t *)graph.xadj().data(),
             (idx_t *)graph.adjncy().data(), (idx_t *)graph.vwgt().data(),
             (idx_t *)0, 0, (idx_t *)&numPart, 0, 0,
             (idx_t *)metisOptions, (idx_t *)&objval, (idx_t *)epart.data());
        }
    
        switch(metisError){
        case METIS_OK:
          break;
        case METIS_ERROR_INPUT:
          Msg::Error("METIS input error");
          return 1;
        case METIS_ERROR_MEMORY:
          Msg::Error("METIS memoty error");
          return 1;
        case METIS_ERROR:
        default:
          Msg::Error("METIS error");
          return 1;
        }
    
        // Check and correct the topology
        for(unsigned int i = 0; i < graph.ne(); i++){
          if(graph.element(i)->getDim() == (int)graph.dim()) continue;
    
          for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
            if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-1){
              if(epart[i] != epart[graph.adjncy(j)]){
                epart[i] = epart[graph.adjncy(j)];
                break;
              }
            }
          }
    
          for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
            if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-2){
              if(epart[i] != epart[graph.adjncy(j)]){
                epart[i] = epart[graph.adjncy(j)];
                break;
              }
            }
          }
    
          for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
            if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-3){
              if(epart[i] != epart[graph.adjncy(j)]){
                epart[i] = epart[graph.adjncy(j)];
                break;
              }
            }
          }
        }
        graph.partition(epart);
    
        Msg::Info("%d partitions, %d total edge-cuts", numPart, objval);
      }
      catch(...) {
        Msg::Error("METIS exception");
        return 2;
      }
    #endif
    
      return 0;
    }
    
    template <class ITERATOR>
    static void assignElementsToEntities(GModel *const model,
                                         hashmap<MElement*, unsigned int> &elmToPartition,
                                         std::vector<partitionRegion *> &newRegions,
                                         ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        const unsigned int partition = elmToPartition[(*it)] - 1;
    
        if(!newRegions[partition]){
          std::vector<unsigned int> partitions;
          partitions.push_back(partition+1);
          partitionRegion *dr = new partitionRegion(model, ++elementaryNumber, partitions);
          model->add(dr);
          newRegions[partition] = dr;
        }
    
        newRegions[partition]->addElement((*it)->getType(), *it);
      }
    }
    
    template <class ITERATOR>
    static void assignElementsToEntities(GModel *const model,
                                         hashmap<MElement*, unsigned int> &elmToPartition,
                                         std::vector<partitionFace *> &newFaces,
                                         ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        const unsigned int partition = elmToPartition[(*it)] - 1;
    
        if(!newFaces[partition]){
          std::vector<unsigned int> partitions;
          partitions.push_back(partition+1);
          partitionFace *df = new partitionFace(model, ++elementaryNumber, partitions);
          model->add(df);
          newFaces[partition] = df;
        }
    
        newFaces[partition]->addElement((*it)->getType(), *it);
      }
    }
    
    template <class ITERATOR>
    static void assignElementsToEntities(GModel *const model,
                                         hashmap<MElement*, unsigned int> &elmToPartition,
                                         std::vector<partitionEdge *> &newEdges,
                                         ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        const unsigned int partition = elmToPartition[(*it)] - 1;
    
        if(!newEdges[partition]){
          std::vector<unsigned int> partitions;
          partitions.push_back(partition+1);
          partitionEdge *de = new partitionEdge(model, ++elementaryNumber, 0, 0, partitions);
          model->add(de);
          newEdges[partition] = de;
        }
    
        newEdges[partition]->addElement((*it)->getType(), *it);
      }
    }
    
    template <class ITERATOR>
    static void assignElementsToEntities(GModel *const model,
                                         hashmap<MElement*, unsigned int> &elmToPartition,
                                         std::vector<partitionVertex *> &newVertices,
                                         ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        const unsigned int partition = elmToPartition[(*it)] - 1;
    
        if(!newVertices[partition]){
          std::vector<unsigned int> partitions;
          partitions.push_back(partition+1);
          partitionVertex *dv = new partitionVertex(model, ++elementaryNumber, partitions);
          model->add(dv);
          newVertices[partition] = dv;
        }
    
        newVertices[partition]->addElement((*it)->getType(), *it);
      }
    }
    
    template <class ITERATOR>
    void setVerticesToEntity(GEntity *const entity, ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
          if(!(*it)->getVertex(i)->onWhat()){
            (*it)->getVertex(i)->setEntity(entity);
            entity->addMeshVertex((*it)->getVertex(i));
          }
        }
      }
    }
    
    // Assign the vertices to its corresponding entity
    static void AssignMeshVertices(GModel *model)
    {
      // Loop over vertices
      for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
        for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
          for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
            (*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
          }
        }
        (*it)->mesh_vertices.clear();
      }
    
      // Loop over edges
      for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
        for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
          for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
            (*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
          }
        }
        (*it)->mesh_vertices.clear();
      }
    
      // Loop over faces
      for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
        for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
          for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
            (*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
          }
        }
        (*it)->mesh_vertices.clear();
      }
    
      // Loop over regions
      for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
        for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
          for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
            (*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
          }
        }
        (*it)->mesh_vertices.clear();
      }
    
    
      // Loop over vertices
      for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
        setVerticesToEntity(*it, (*it)->points.begin(), (*it)->points.end());
      }
    
      // Loop over edges
      for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
          setVerticesToEntity(*it, (*it)->lines.begin(), (*it)->lines.end());
      }
    
      // Loop over faces
      for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
          setVerticesToEntity(*it, (*it)->triangles.begin(), (*it)->triangles.end());
          setVerticesToEntity(*it, (*it)->quadrangles.begin(), (*it)->quadrangles.end());
      }
    
      // Loop over regions
      for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
          setVerticesToEntity(*it, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
          setVerticesToEntity(*it, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
          setVerticesToEntity(*it, (*it)->prisms.begin(), (*it)->prisms.end());
          setVerticesToEntity(*it, (*it)->pyramids.begin(), (*it)->pyramids.end());
          setVerticesToEntity(*it, (*it)->trihedra.begin(), (*it)->trihedra.end());
      }
    }
    
    static void fillConnectedElements(std::vector< std::set<MElement*> > &connectedElements,
                                      Graph &graph)
    {
      std::stack<unsigned int> elementStack;
      std::set<MElement*> elements;
      unsigned int startElement = 0;
      bool stop = true;
      unsigned int size = 0;
      int isolatedElements = 0;
    
      do {
        //Inititalization
        elementStack.push(startElement);
        elements.insert(graph.element(startElement));
    
        while(elementStack.size() != 0){
          unsigned int top = elementStack.top();
          elementStack.pop();
          elements.insert(graph.element(top));
    
          for(unsigned int i = graph.xadj(top); i < graph.xadj(top+1); i++){
            if(graph.adjncy(i) != 0){
              elementStack.push(graph.adjncy(i));
              graph.adjncy(i,0);
            }
          }
        }
        connectedElements.push_back(elements);
        size += elements.size();
        elements.clear();
    
        stop = (size == graph.ne() ? true : false);
    
        startElement = 0;
        if(!stop){
          for(unsigned int i = 0; i < graph.ne(); i++){
            for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
              if(graph.adjncy(j) != 0){
                startElement = i;
                i = graph.ne();
                break;
              }
            }
          }
          if(startElement == 0){
            int skipIsolatedElements = 0;
            for(unsigned int i = 1; i < graph.ne(); i++){
              if(graph.xadj(i) == graph.xadj(i+1)){
                if(skipIsolatedElements == isolatedElements){
                  startElement = i;
                  isolatedElements++;
                  break;
                }
                skipIsolatedElements++;
              }
            }
          }
        }
        //break;
      } while(!stop);
    }
    
    static bool dividedNonConnectedEntities(GModel *const model, int dim,
                                            std::set<GRegion*, GEntityLessThan> &regions,
                                            std::set<GFace*, GEntityLessThan> &faces,
                                            std::set<GEdge*, GEntityLessThan> &edges,
                                            std::set<GVertex*, GEntityLessThan> &vertices)
    {
      bool ret = false;
      // Loop over vertices
      if(dim < 0 || dim == 0){
        int elementaryNumber = model->getMaxElementaryNumber(0);
    
        for(GModel::const_viter it = vertices.begin(); it != vertices.end(); ++it){
          if((*it)->geomType() == GEntity::PartitionVertex){
            partitionVertex* vertex = static_cast<partitionVertex*>(*it);
    
            if(vertex->getNumMeshElements() > 1){
              ret = true;
              for(unsigned int i = 0; i < vertex->getNumMeshElements(); i++){
                // Create the new partitionVertex
                partitionVertex *pvertex = new partitionVertex
                (model, ++elementaryNumber, vertex->getPartitions());
                // Assign physicals and parent entity
                std::vector<int> physicalTags = vertex->getPhysicalEntities();
                for(unsigned int j = 0; j < physicalTags.size(); j++){
                  pvertex->addPhysicalEntity(physicalTags[j]);
                }
                pvertex->setParentEntity(vertex->getParentEntity());
                // Add to model
                model->add(pvertex);
                // Add elements
                pvertex->addElement(vertex->getMeshElement(i)->getType(),
                                    vertex->getMeshElement(i));
                // Move B-Rep
                std::vector<GEdge*> BRepEdges = vertex->edges();
                if(!BRepEdges.empty()){
                  for(std::vector<GEdge*>::iterator itBRep = BRepEdges.begin();
                      itBRep != BRepEdges.end(); ++itBRep){
                    if(vertex == (*itBRep)->getBeginVertex()){
                      (*itBRep)->setVertex(pvertex, 1);
                      pvertex->addEdge(*itBRep);
                    }
                    if(vertex == (*itBRep)->getEndVertex()){
                      (*itBRep)->setVertex(pvertex, -1);
                      pvertex->addEdge(*itBRep);
                    }
                  }
                }
              }
    
              model->remove(vertex);
              vertex->points.clear();
              vertex->mesh_vertices.clear();
              delete vertex;
            }
          }
        }
      }
    
      // Loop over edges
      if(dim < 0 || dim == 1){
        // We build a graph
        Graph graph(model);
        graph.ne(model->getNumMeshElements(1));
        graph.nn(model->getNumMeshVertices(1));
        graph.dim(model->getDim());
        graph.elementResize(graph.ne());
        graph.vertexResize(model->getMaxVertexNumber());
        graph.eptrResize(graph.ne()+1);
        graph.eptr(0,0);
        const int eindSize = getSizeOfEind(model);
        graph.eindResize(eindSize);
    
        int elementaryNumber = model->getMaxElementaryNumber(1);
    
        for(GModel::const_eiter it = edges.begin(); it != edges.end(); ++it){
          if((*it)->geomType() == GEntity::PartitionCurve){
            partitionEdge* edge = static_cast<partitionEdge*>(*it);
    
            graph.ne(edge->getNumMeshElements());
            graph.dim(1);
            graph.eptr(0,0);
            graph.clearDualGraph();
            graph.eraseVertex();
    
            size_t eptrIndex = 0;
            size_t eindIndex = 0;
            size_t numVertex = 0;
    
            fillElementsToNodesMap(graph, edge, eptrIndex, eindIndex, numVertex,
                                   edge->lines.begin(), edge->lines.end());
            graph.nn(numVertex);
            graph.createDualGraph(false);
    
            // if a graph contains more than ((n-1)*(n-2))/2 edges
            // (where n is the number of nodes), then it is connected.
            if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
              continue;
            }
    
            std::vector< std::set<MElement*> > connectedElements;
            fillConnectedElements(connectedElements, graph);
    
            if(connectedElements.size() > 1){
              ret = true;
              std::vector<GFace*> BRepFaces = edge->faces();
    
              std::vector<int> oldOrientations;
              oldOrientations.reserve(BRepFaces.size());
    
              if(!BRepFaces.empty()){
                for(std::vector<GFace*>::iterator itBRep = BRepFaces.begin();
                    itBRep !=  BRepFaces.end(); ++itBRep){
                  oldOrientations.push_back((*itBRep)->delEdge(edge));
                }
              }
    
              for(unsigned int i = 0; i < connectedElements.size(); i++){
                // Create the new partitionEdge
                partitionEdge *pedge = new partitionEdge
                (model, ++elementaryNumber, 0, 0, edge->getPartitions());
                // Assign physicals and parent entity
                std::vector<int> physicalTags = edge->getPhysicalEntities();
                for(unsigned int j = 0; j < physicalTags.size(); j++){
                  pedge->addPhysicalEntity(physicalTags[j]);
                }
                pedge->setParentEntity(edge->getParentEntity());
                // Add to model
                model->add(pedge);
                for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
                    itSet != connectedElements[i].end(); ++itSet){
                  // Add elements
                  pedge->addElement((*itSet)->getType(), (*itSet));
                }
                // Move B-Rep
                if(BRepFaces.size() > 0){
                  int i = 0;
                  for(std::vector<GFace*>::iterator itBRep = BRepFaces.begin();
                      itBRep !=  BRepFaces.end(); ++itBRep){
                    (*itBRep)->setEdge(pedge, oldOrientations[i]);
                    pedge->addFace(*itBRep);
                    i++;
                  }
                }
              }
    
              model->remove(edge);
              edge->lines.clear();
              edge->mesh_vertices.clear();
              delete edge;
            }
    
            connectedElements.clear();
          }
        }
      }
    
      // Loop over faces
      if(dim < 0 || dim == 2){
        // We build a graph
        Graph graph(model);
        graph.ne(model->getNumMeshElements(2));
        graph.nn(model->getNumMeshVertices(2));
        graph.dim(model->getDim());
        graph.elementResize(graph.ne());
        graph.vertexResize(model->getMaxVertexNumber());
        graph.eptrResize(graph.ne()+1);
        graph.eptr(0,0);
        const int eindSize = getSizeOfEind(model);
        graph.eindResize(eindSize);
    
        int elementaryNumber = model->getMaxElementaryNumber(2);
    
        for(GModel::const_fiter it = faces.begin(); it != faces.end(); ++it){
          if((*it)->geomType() == GEntity::PartitionSurface){
            partitionFace* face = static_cast<partitionFace*>(*it);
    
            graph.ne(face->getNumMeshElements());
            graph.dim(2);
            graph.eptr(0,0);
            graph.clearDualGraph();
            graph.eraseVertex();
    
            size_t eptrIndex = 0;
            size_t eindIndex = 0;
            size_t numVertex = 0;
    
            fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
                                   face->triangles.begin(), face->triangles.end());
            fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
                                   face->quadrangles.begin(), face->quadrangles.end());
            graph.nn(numVertex);
            graph.createDualGraph(false);
    
            // if a graph contains more than ((n-1)*(n-2))/2 edges
            // (where n is the number of nodes), then it is connected.
            if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
              continue;
            }
    
            std::vector< std::set<MElement*> > connectedElements;
            fillConnectedElements(connectedElements, graph);
    
            if(connectedElements.size() > 1){
              ret = true;
              std::list<GRegion*> BRepRegions = face->regions();
              std::vector<int> oldOrientations;
              if(BRepRegions.size() > 0){
                for(std::list<GRegion*>::iterator itBRep = BRepRegions.begin();
                    itBRep !=  BRepRegions.end(); ++itBRep){
                  oldOrientations.push_back((*itBRep)->delFace(face));
                }
              }
    
              for(unsigned int i = 0; i < connectedElements.size(); i++){
                // Create the new partitionFace
                partitionFace *pface = new partitionFace
                (model, ++elementaryNumber, face->getPartitions());
                // Assign physicals and parent entity
                std::vector<int> physicalTags = face->getPhysicalEntities();
                for(unsigned int j = 0; j < physicalTags.size(); j++){
                  pface->addPhysicalEntity(physicalTags[j]);
                }
                pface->setParentEntity(face->getParentEntity());
                // Add to model
                model->add(pface);
                for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
                    itSet != connectedElements[i].end(); ++itSet){
                  // Add elements
                  pface->addElement((*itSet)->getType(), (*itSet));
                }
                // Move B-Rep
                if(BRepRegions.size() > 0){
                  int i = 0;
                  for(std::list<GRegion*>::iterator itBRep = BRepRegions.begin();
                      itBRep !=  BRepRegions.end(); ++itBRep){
                    (*itBRep)->setFace(pface, oldOrientations[i]);
                    pface->addRegion(*itBRep);
                    i++;
                  }
                }
              }
    
              model->remove(face);
              face->triangles.clear();
              face->quadrangles.clear();
              face->mesh_vertices.clear();
              delete face;
            }
    
            connectedElements.clear();
          }
        }
      }
    
      // Loop over regions
      if(dim < 0 || dim == 3){
        // We build a graph
        Graph graph(model);
        graph.ne(model->getNumMeshElements(3));
        graph.nn(model->getNumMeshVertices(3));
        graph.dim(model->getDim());
        graph.elementResize(graph.ne());
        graph.vertexResize(model->getMaxVertexNumber());
        graph.eptrResize(graph.ne()+1);
        graph.eptr(0,0);
        const int eindSize = getSizeOfEind(model);
        graph.eindResize(eindSize);
    
        int elementaryNumber = model->getMaxElementaryNumber(3);
    
        for(GModel::const_riter it = regions.begin(); it != regions.end(); ++it){
          if((*it)->geomType() == GEntity::PartitionVolume){
            partitionRegion* region = static_cast<partitionRegion*>(*it);
    
            graph.ne(region->getNumMeshElements());
            graph.dim(3);
            graph.eptr(0,0);
            graph.clearDualGraph();
            graph.eraseVertex();
    
            size_t eptrIndex = 0;
            size_t eindIndex = 0;
            size_t numVertex = 0;
    
            fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
                                   region->tetrahedra.begin(), region->tetrahedra.end());
            fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
                                   region->hexahedra.begin(), region->hexahedra.end());
            fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
                                   region->prisms.begin(), region->prisms.end());
            fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
                                   region->pyramids.begin(), region->pyramids.end());
            fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
                                   region->trihedra.begin(), region->trihedra.end());
            graph.nn(numVertex);
            graph.createDualGraph(false);
    
            // if a graph contains more than ((n-1)*(n-2))/2 edges
            // (where n is the number of nodes), then it is connected.
            if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
              continue;
            }
    
            std::vector< std::set<MElement*> > connectedElements;
            fillConnectedElements(connectedElements, graph);
    
            if(connectedElements.size() > 1){
              ret = true;
              for(unsigned int i = 0; i < connectedElements.size(); i++){
                // Create the new partitionRegion
                partitionRegion *pregion = new partitionRegion
                (model, ++elementaryNumber, region->getPartitions());
                // Assign physicals and parent entity
                std::vector<int> physicalTags = (*it)->getPhysicalEntities();
                for(unsigned int j = 0; j < physicalTags.size(); j++){
                  pregion->addPhysicalEntity(physicalTags[j]);
                }
                pregion->setParentEntity(region->getParentEntity());
                // Add to model
                model->add(pregion);
                for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
                    itSet != connectedElements[i].end(); ++itSet){
                  //Add elements
                  pregion->addElement((*itSet)->getType(), (*itSet));
                }
              }
    
              model->remove(region);
              region->tetrahedra.clear();
              region->hexahedra.clear();
              region->prisms.clear();
              region->pyramids.clear();
              region->trihedra.clear();
              region->mesh_vertices.clear();
              delete region;
            }
    
            connectedElements.clear();
          }
        }
      }
    
      return ret;
    }
    
    // Create the new volume entities (omega)
    static void CreateNewEntities(GModel *const model,
                                  hashmap<MElement*, unsigned int> &elmToPartition)
    {
      std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
      std::set<GFace*, GEntityLessThan> faces = model->getFaces();
      std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
      std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
    
      int elementaryNumber = model->getMaxElementaryNumber(0);
      for(GModel::const_viter it = vertices.begin(); it != vertices.end(); ++it){
        std::vector<partitionVertex *> newVertices(model->getNumPartitions(), 0);
    
        assignElementsToEntities(model, elmToPartition, newVertices,
                                 (*it)->points.begin(), (*it)->points.end(), elementaryNumber);
    
        for(unsigned int i = 0; i < model->getNumPartitions(); i++){
          if(newVertices[i]){
            static_cast<partitionVertex*>(newVertices[i])->setParentEntity((*it));
            std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
            for(unsigned int j = 0; j < physicalEntities.size(); j++){
              newVertices[i]->addPhysicalEntity(physicalEntities[j]);
            }
          }
        }
    
        (*it)->mesh_vertices.clear();
    
        (*it)->points.clear();
      }
    
      elementaryNumber = model->getMaxElementaryNumber(1);
      for(GModel::const_eiter it = edges.begin(); it != edges.end(); ++it){
        std::vector<partitionEdge *> newEdges(model->getNumPartitions(), 0);
    
        assignElementsToEntities(model, elmToPartition, newEdges,
                                 (*it)->lines.begin(), (*it)->lines.end(), elementaryNumber);
    
        for(unsigned int i = 0; i < model->getNumPartitions(); i++){
          if(newEdges[i]){
            static_cast<partitionEdge*>(newEdges[i])->setParentEntity(*it);
            std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
            for(unsigned int j = 0; j < physicalEntities.size(); j++){
              newEdges[i]->addPhysicalEntity(physicalEntities[j]);
            }
          }
        }
    
        (*it)->mesh_vertices.clear();
    
        (*it)->lines.clear();
      }
    
      elementaryNumber = model->getMaxElementaryNumber(2);
      for(GModel::const_fiter it = faces.begin(); it != faces.end(); ++it){
        std::vector<partitionFace *> newFaces(model->getNumPartitions(), 0);
    
        assignElementsToEntities
        (model, elmToPartition, newFaces,
         (*it)->triangles.begin(), (*it)->triangles.end(), elementaryNumber);
        assignElementsToEntities
        (model, elmToPartition, newFaces,
         (*it)->quadrangles.begin(), (*it)->quadrangles.end(), elementaryNumber);
    
        std::list<GRegion*> BRepRegions = (*it)->regions();
        for(unsigned int i = 0; i < model->getNumPartitions(); i++){
          if(newFaces[i]){
            static_cast<partitionFace*>(newFaces[i])->setParentEntity(*it);
            std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
            for(unsigned int j = 0; j < physicalEntities.size(); j++){
              newFaces[i]->addPhysicalEntity(physicalEntities[j]);
            }
          }
        }
    
        (*it)->mesh_vertices.clear();
    
        (*it)->triangles.clear();
        (*it)->quadrangles.clear();
      }
    
      elementaryNumber = model->getMaxElementaryNumber(3);
      for(GModel::const_riter it = regions.begin(); it != regions.end(); ++it){
        std::vector<partitionRegion *> newRegions(model->getNumPartitions(), 0);
    
        assignElementsToEntities
          (model, elmToPartition, newRegions,
           (*it)->tetrahedra.begin(), (*it)->tetrahedra.end(), elementaryNumber);
        assignElementsToEntities
          (model, elmToPartition, newRegions,
           (*it)->hexahedra.begin(), (*it)->hexahedra.end(), elementaryNumber);
        assignElementsToEntities
          (model, elmToPartition, newRegions,
           (*it)->prisms.begin(), (*it)->prisms.end(), elementaryNumber);
        assignElementsToEntities
          (model, elmToPartition, newRegions,
           (*it)->pyramids.begin(), (*it)->pyramids.end(), elementaryNumber);
        assignElementsToEntities
          (model, elmToPartition, newRegions,
           (*it)->trihedra.begin(), (*it)->trihedra.end(), elementaryNumber);
    
        for(unsigned int i = 0; i < model->getNumPartitions(); i++){
          if(newRegions[i]){
            static_cast<partitionRegion*>(newRegions[i])->setParentEntity(*it);
            std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
            for(unsigned int j = 0; j < physicalEntities.size(); j++){
              newRegions[i]->addPhysicalEntity(physicalEntities[j]);
            }
          }
        }
    
        (*it)->mesh_vertices.clear();
    
        (*it)->tetrahedra.clear();
        (*it)->hexahedra.clear();
        (*it)->prisms.clear();
        (*it)->pyramids.clear();
        (*it)->trihedra.clear();
      }
    
      // If we don't create the partition topology let's just assume that the user
      // does not care about multi-connected partitions or partition boundaries.
      if(!CTX::instance()->mesh.partitionCreateTopology) return;
      regions = model->getRegions();
      faces = model->getFaces();
      edges = model->getEdges();
      vertices = model->getVertices();
      dividedNonConnectedEntities(model, -1, regions, faces, edges, vertices);
    }
    
    static void fillElementToEntity(GModel *const model,
                                    hashmap<MElement*, GEntity*> &elmToEntity, int dim)
    {
      // Loop over regions
      if(dim < 0 || dim == 3){
        for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
          for(std::vector<MTetrahedron*>::iterator itElm = (*it)->tetrahedra.begin();
              itElm != (*it)->tetrahedra.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
          for(std::vector<MHexahedron*>::iterator itElm = (*it)->hexahedra.begin();
              itElm != (*it)->hexahedra.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
          for(std::vector<MPrism*>::iterator itElm = (*it)->prisms.begin();
              itElm != (*it)->prisms.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
          for(std::vector<MPyramid*>::iterator itElm = (*it)->pyramids.begin();
              itElm != (*it)->pyramids.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
          for(std::vector<MTrihedron*>::iterator itElm = (*it)->trihedra.begin();
              itElm != (*it)->trihedra.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
        }
      }
    
      // Loop over faces
      if(dim < 0 || dim == 2){
        for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
          for(std::vector<MTriangle*>::iterator itElm = (*it)->triangles.begin();
              itElm != (*it)->triangles.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
          for(std::vector<MQuadrangle*>::iterator itElm = (*it)->quadrangles.begin();
              itElm != (*it)->quadrangles.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
        }
      }
    
      // Loop over edges
      if(dim < 0 || dim == 1){
        for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
          for(std::vector<MLine*>::iterator itElm = (*it)->lines.begin();
              itElm != (*it)->lines.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
        }
      }
    
      // Loop over vertices
      if(dim < 0 || dim == 0){
        for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
          for(std::vector<MPoint*>::iterator itElm = (*it)->points.begin();
              itElm != (*it)->points.end(); ++itElm)
            elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
        }
      }
    }
    
    static MElement* getReferenceElement(const std::vector< std::pair<MElement*,
                                         std::vector<unsigned int> > > &elementPairs)
    {
      unsigned int min = std::numeric_limits<unsigned int>::max();
      std::vector< std::pair<MElement*, std::vector<unsigned int> > > minSizeElementPairs;
      std::vector< std::pair<MElement*, std::vector<unsigned int> > > minSizeElementPairsTmp;
    
      // Take only the elements having the less partition in commun. For exemple we
      // take (1,2) and (3,8) but not (2,5,9) or (1,4,5,7)
      for(unsigned int i = 0; i < elementPairs.size(); i++){
        if(min > elementPairs[i].second.size()){
          minSizeElementPairs.clear();
          min = elementPairs[i].second.size();
          minSizeElementPairs.push_back(elementPairs[i]);
        }
        else if(min == elementPairs[i].second.size()){
          minSizeElementPairs.push_back(elementPairs[i]);
        }
      }
    
      // Check if the element separated different partitions
      if(minSizeElementPairs.size() == elementPairs.size()){
        bool isEqual = true;
        for(unsigned int i = 1; i < minSizeElementPairs.size(); i++){
          if(minSizeElementPairs[i].second != minSizeElementPairs[0].second){
            isEqual = false;
            break;
          }
        }
        if(isEqual) return 0;
      }
    
      while(minSizeElementPairs.size() > 1){
        min = std::numeric_limits<unsigned int>::max();
        for(unsigned int i = 0; i < minSizeElementPairs.size(); i++){
          // The partition vector is sorted thus we can check only the first element
          if(minSizeElementPairs[i].second.size() == 0) return minSizeElementPairs[0].first;
          if(min > minSizeElementPairs[i].second[0]){
            min = minSizeElementPairs[i].second[0];
          }
        }
    
        for(unsigned int i = 0; i < minSizeElementPairs.size(); i++){
          if(min == minSizeElementPairs[i].second[0]){
            minSizeElementPairs[i].second.erase(minSizeElementPairs[i].second.begin());
            minSizeElementPairsTmp.push_back(minSizeElementPairs[i]);
          }
        }
        minSizeElementPairs.clear();
    #if __cplusplus >= 201103L
        minSizeElementPairs = std::move(minSizeElementPairsTmp);
    #else
        minSizeElementPairs = minSizeElementPairsTmp;
    #endif
        minSizeElementPairsTmp.clear();
      }
    
      return minSizeElementPairs[0].first;
    }
    
    static void getPartitionInVector(std::vector<unsigned int> &partitions,
                                     const std::vector< std::pair<MElement*,
                                     std::vector<unsigned int> > > &boundaryPair)
    {
      for(unsigned int i = 0; i < boundaryPair.size(); i++){
        for(unsigned int j = 0; j < boundaryPair[i].second.size(); j++){
          if(std::find(partitions.begin(), partitions.end(), boundaryPair[i].second[j]) ==
             partitions.end()){
            partitions.push_back(boundaryPair[i].second[j]);
          }
        }
      }
    
      std::sort(partitions.begin(), partitions.end());
    }
    
    static partitionFace *assignPartitionBoundary(
      GModel *const model, MFace &me, MElement *reference,
      const std::vector<unsigned int> &partitions,
      std::multimap<partitionFace *, GEntity *, Less_partitionFace> &pfaces,
      hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
    {
      partitionFace *newEntity = 0;
      partitionFace pf(model, 1, partitions);
      std::pair< std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator,
                 std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator>
        ret = pfaces.equal_range(&pf);
    
      partitionFace *ppf = 0;
      // Create the new partition entity for the mesh
      if(ret.first == ret.second){
        // Create new entity and add them to the model
        ppf = new partitionFace(model, ++numEntity, partitions);
        pfaces.insert(std::pair<partitionFace*, GEntity*>(ppf, elementToEntity[reference]));
        model->add(ppf);
        newEntity = ppf;
      }
      else{
        for(std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator it = ret.first;
            it != ret.second; ++it){
          if(elementToEntity[reference] == (*it).second){
            ppf = (*it).first;
          }
        }
        if(!ppf){
          // Create new entity and add them to the model
          ppf = new partitionFace(model, ++numEntity, partitions);
          pfaces.insert(std::pair<partitionFace*, GEntity*>(ppf, elementToEntity[reference]));
          model->add(ppf);
          newEntity = ppf;
        }
      }
    
      int numFace = 0;
      for(int i = 0; i < reference->getNumFaces(); i++){
        if(reference->getFace(i) == me){
          numFace = i;
          break;
        }
      }
    
      if(me.getNumVertices() == 3){
        std::vector<MVertex*> verts;
        reference->getFaceVertices(numFace, verts);
    
        if(verts.size() == 3){
          MTriangle *element = new MTriangle(verts);
          ppf->addTriangle(element);
        }
        else if(verts.size() == 6){
          MTriangle6 *element = new MTriangle6(verts);
          ppf->addTriangle(element);
        }
        else{
          MTriangleN *element = new MTriangleN(verts, verts[0]->getPolynomialOrder());
          ppf->addTriangle(element);
        }
      }
      else if(me.getNumVertices() == 4){
        std::vector<MVertex*> verts;
        reference->getFaceVertices(numFace, verts);
    
        if(verts.size() == 4){
          MQuadrangle *element = new MQuadrangle(verts);
          ppf->addQuadrangle(element);
        }
        else if(verts.size() == 8){
          MQuadrangle8 *element = new MQuadrangle8(verts);
          ppf->addQuadrangle(element);
        }
        else if(verts.size() == 9){
          MQuadrangle9 *element = new MQuadrangle9(verts);
          ppf->addQuadrangle(element);
        }
        else{
          MQuadrangleN *element = new MQuadrangleN(verts, verts[0]->getPolynomialOrder());
          ppf->addQuadrangle(element);
        }
      }
    
      return newEntity;
    }
    
    static partitionEdge *assignPartitionBoundary(
      GModel *const model, MEdge &me, MElement *reference,
      const std::vector<unsigned int> &partitions,
      std::multimap<partitionEdge *, GEntity *, Less_partitionEdge> &pedges,
      hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
    {
      partitionEdge *newEntity = 0;
      partitionEdge pe(model, 1, 0, 0, partitions);
      std::pair< std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator,
                 std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator>
        ret = pedges.equal_range(&pe);
    
      partitionEdge *ppe = 0;
      // Create the new partition entity for the mesh
      if(ret.first == ret.second){
        // Create new entity and add them to the model
        ppe = new partitionEdge(model, ++numEntity, 0, 0, partitions);
        pedges.insert(std::pair<partitionEdge*, GEntity*>(ppe, elementToEntity[reference]));
        model->add(ppe);
        newEntity = ppe;
      }
      else{
        for(std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator it = ret.first;
            it != ret.second; ++it){
          if(elementToEntity[reference] == (*it).second){
            ppe = (*it).first;
          }
        }
        if(!ppe){
          // Create new entity and add them to the model
          ppe = new partitionEdge(model, ++numEntity, 0, 0, partitions);
          pedges.insert(std::pair<partitionEdge*, GEntity*>(ppe, elementToEntity[reference]));
          model->add(ppe);
          newEntity = ppe;
        }
      }
    
      int numEdge = 0;
      for(int i = 0; i < reference->getNumEdges(); i++){
        if(reference->getEdge(i) == me){
          numEdge = i;
          break;
        }
      }
    
      if(me.getNumVertices() == 2){
        std::vector<MVertex*> verts;
        reference->getEdgeVertices(numEdge, verts);
    
        if(verts.size() == 2){
          MLine *element = new MLine(verts);
          ppe->addLine(element);
        }
        else if(verts.size() == 3){
          MLine3 *element = new MLine3(verts);
          ppe->addLine(element);
        }
        else{
          MLineN *element = new MLineN(verts);
          ppe->addLine(element);
        }
      }
    
      return newEntity;
    }
    
    static partitionVertex *assignPartitionBoundary(
      GModel *const model, MVertex *ve, MElement *reference,
      const std::vector<unsigned int> &partitions,
      std::multimap<partitionVertex *, GEntity *, Less_partitionVertex> &pvertices,
      hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
    {
      partitionVertex *newEntity = 0;
      partitionVertex pv(model, 1, partitions);
      std::pair< std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator,
                 std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator >
        ret = pvertices.equal_range(&pv);
    
      partitionVertex *ppv = 0;
      // Create the new partition entity for the mesh
      if(ret.first == ret.second){
        ppv = new partitionVertex(model, ++numEntity, partitions);
        pvertices.insert(std::pair<partitionVertex*, GEntity*>(ppv, elementToEntity[reference]));
        model->add(ppv);
        newEntity = ppv;
      }
      else{
        for(std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator it =
              ret.first; it != ret.second; ++it){
          if(elementToEntity[reference] == (*it).second){
            ppv = (*it).first;
          }
        }
        if(!ppv){
          // Create new entity and add them to the model
          ppv = new partitionVertex(model, ++numEntity, partitions);
          pvertices.insert(std::pair<partitionVertex*, GEntity*>
                           (ppv, elementToEntity[reference]));
          model->add(ppv);
          newEntity = ppv;
        }
      }
    
      ppv->addPoint(new MPoint(ve));
    
      return newEntity;
    }
    
    static int computeOrientation(MElement* reference, MElement* element)
    {
      if(element->getDim() == 2)
      {
        std::vector<MVertex*> vertices;
        element->getVertices(vertices);
        MFace face = element->getFace(0);
        for(int i = 0; i < reference->getNumFaces(); i++){
          if(reference->getFace(i) == face){
            std::vector<MVertex*> referenceVertices;
            reference->getFaceVertices(i, referenceVertices);
    
            if(referenceVertices == vertices) return 1;
            else return -1;
          }
        }
      }
      else if(element->getDim() == 1)
      {
        std::vector<MVertex*> vertices;
        element->getVertices(vertices);
        MEdge face = element->getEdge(0);
        for(int i = 0; i < reference->getNumEdges(); i++){
          if(reference->getEdge(i) == face){
            std::vector<MVertex*> referenceVertices;
            reference->getEdgeVertices(i, referenceVertices);
    
            if(referenceVertices == vertices) return 1;
            else return -1;
          }
        }
      }
      else if(element->getDim() == 0)
      {
        std::vector<MVertex*> vertices;
        element->getVertices(vertices);
    
        std::vector<MVertex*> referenceVertices;
        reference->getVertices(referenceVertices);
    
        if(referenceVertices[0] == vertices[0]) return 1;
        if(referenceVertices[1] == vertices[0]) return -1;
      }
    
      return 0;
    }
    
    static void assignBrep(GModel *const model, std::map<GEntity*, MElement*>
                           &boundaryEntityAndRefElement, GEntity *e)
    {
      if(e->dim() == 2){
        partitionFace* entity = static_cast<partitionFace*>(e);
    
        for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
            it != boundaryEntityAndRefElement.end(); ++it){
          static_cast<GRegion*>(it->first)->setFace
            (entity, computeOrientation(it->second, entity->getMeshElement(0)));
          entity->addRegion(static_cast<GRegion*>(it->first));
        }
      }
      else if(e->dim() == 1){
        partitionEdge* entity = static_cast<partitionEdge*>(e);
    
        for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
            it != boundaryEntityAndRefElement.end(); ++it){
          static_cast<GFace*>(it->first)->setEdge
            (entity, computeOrientation(it->second, entity->getMeshElement(0)));
          entity->addFace(static_cast<GFace*>(it->first));
        }
      }
      else if(e->dim() == 0){
        partitionVertex* entity = static_cast<partitionVertex*>(e);
    
        for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
            it != boundaryEntityAndRefElement.end(); ++it){
          static_cast<GEdge*>(it->first)->setVertex
            (entity, computeOrientation(it->second, entity->getMeshElement(0)));
          entity->addEdge(static_cast<GEdge*>(it->first));
        }
      }
    }
    
    void assignNewEntityBRep(Graph &graph, hashmap<MElement*, GEntity*> &elementToEntity)
    {
      std::set< std::pair<GEntity*, GEntity*> > brepWithoutOri;
      hashmap<GEntity*, std::set<std::pair<int, GEntity*> > > brep;
      for(unsigned int i = 0; i < graph.ne(); i++){
        MElement *current = graph.element(i);
        for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++)
        {
          if(current->getDim() == graph.element(graph.adjncy(j))->getDim()+1){
            GEntity *g1 = elementToEntity[current];
            GEntity *g2 = elementToEntity[graph.element(graph.adjncy(j))];
            if(brepWithoutOri.find(std::pair<GEntity*, GEntity*>(g1,g2)) == brepWithoutOri.end()){
              const int ori = computeOrientation(current,graph.element(graph.adjncy(j)));
              brepWithoutOri.insert(std::pair<GEntity*, GEntity*>(g1,g2));
              brep[g1].insert(std::pair<int, GEntity*>(ori,g2));
            }
          }
        }
      }
    
      for(hashmap<GEntity*, std::set<std::pair<int, GEntity*> > >::iterator
            it = brep.begin(); it != brep.end(); ++it){
        switch (it->first->dim()) {
            case 3:
            for(std::set<std::pair<int, GEntity*> >::iterator itSet =
                  it->second.begin(); itSet != it->second.end(); ++itSet){
              static_cast<GRegion*>(it->first)->setFace
                (static_cast<GFace*>(itSet->second), itSet->first);
              static_cast<GFace*>(itSet->second)->addRegion
                (static_cast<GRegion*>(it->first));
            }
            break;
            case 2:
            for(std::set<std::pair<int, GEntity*> >::iterator itSet =
                  it->second.begin(); itSet != it->second.end(); ++itSet){
              static_cast<GFace*>(it->first)->setEdge
                (static_cast<GEdge*>(itSet->second), itSet->first);
              static_cast<GEdge*>(itSet->second)->addFace
                (static_cast<GFace*>(it->first));
            }
            break;
            case 1:
            for(std::set<std::pair<int, GEntity*> >::iterator itSet =
                  it->second.begin(); itSet != it->second.end(); ++itSet){
              static_cast<GEdge*>(it->first)->setVertex
                (static_cast<GVertex*>(itSet->second), itSet->first);
              static_cast<GVertex*>(itSet->second)->addEdge
                (static_cast<GEdge*>(it->first));
            }
            break;
          default:
            break;
        }
      }
    }
    
    // Create the new entities between each partitions (sigma and bndSigma).
    static void CreatePartitionTopology(GModel *const model,
                                        const std::vector< std::set<MElement*> >
                                        &boundaryElements, Graph &meshGraph)
    {
      const int meshDim = model->getDim();
      hashmap<MElement*, GEntity*> elementToEntity;
      fillElementToEntity(model, elementToEntity, -1);
      assignNewEntityBRep(meshGraph, elementToEntity);
    
      std::multimap<partitionFace*, GEntity*, Less_partitionFace> pfaces;
      std::multimap<partitionEdge*, GEntity*, Less_partitionEdge> pedges;
      std::multimap<partitionVertex*, GEntity*, Less_partitionVertex> pvertices;
    
      hashmapface faceToElement;
      hashmapedge edgeToElement;
      hashmapvertex vertexToElement;
    
      std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
      std::set<GFace*, GEntityLessThan> faces = model->getFaces();
      std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
      std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
    
      if (meshDim >= 3){ // Create partition faces
        Msg::Info(" - Creating partition faces");
    
        for(unsigned int i = 0; i < model->getNumPartitions(); i++){
          for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
              it != boundaryElements[i].end(); ++it){
            for(int j = 0; j < (*it)->getNumFaces(); j++){
              faceToElement[(*it)->getFace(j)].push_back
              (std::pair<MElement*, std::vector<unsigned int> >
               (*it, std::vector<unsigned int>(1, i+1)));
            }
          }
        }
        int numFaceEntity = model->getMaxElementaryNumber(2);
        for(hashmapface::const_iterator it = faceToElement.begin();
            it != faceToElement.end(); ++it){
          MFace f = it->first;
    
          std::vector<unsigned int> partitions;
          getPartitionInVector(partitions, it->second);
          if(partitions.size() < 2) continue;
    
          MElement* reference = getReferenceElement(it->second);
          if(!reference) continue;
    
          partitionFace *pf = assignPartitionBoundary
            (model, f, reference, partitions, pfaces, elementToEntity, numFaceEntity);
          if(pf){
            std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
            for(unsigned int i = 0; i < it->second.size(); i++)
              boundaryEntityAndRefElement.insert
                (std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
                                                it->second[i].first));
    
            assignBrep(model, boundaryEntityAndRefElement, pf);
          }
        }
        faceToElement.clear();
    
        faces = model->getFaces();
        dividedNonConnectedEntities(model, 2, regions, faces, edges, vertices);
        elementToEntity.clear();
        fillElementToEntity(model, elementToEntity, 2);
      }
    
      if (meshDim >= 2){ // Create partition edges
        Msg::Info(" - Creating partition edges");
    
        if (meshDim == 2){
          for(unsigned int i = 0; i < model->getNumPartitions(); i++){
            for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
                it != boundaryElements[i].end(); ++it){
              for(int j = 0; j < (*it)->getNumEdges(); j++){
                edgeToElement[(*it)->getEdge(j)].push_back
                (std::pair<MElement*, std::vector<unsigned int> >
                 (*it, std::vector<unsigned int>(1, i+1)));
              }
            }
          }
        }
        else{
          Graph subGraph(model);
          MakeGraph(model, subGraph, 2);
          subGraph.createDualGraph(false);
          std::vector<unsigned int>part(subGraph.ne());
          int partIndex = 0;
    
          std::map<unsigned int, std::vector<unsigned int> > mapOfPartitions;
          unsigned int mapOfPartitionsTag = 0;
          for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
            if((*it)->geomType() == GEntity::PartitionSurface){
              std::vector<unsigned int> partitions =
                static_cast<partitionFace*>(*it)->getPartitions();
              mapOfPartitions.insert(std::pair<unsigned int, std::vector<unsigned int> >
                                     (mapOfPartitionsTag, partitions));
              // Must absolutely be in the same order as in the MakeGraph function
              for(std::vector<MTriangle*>::iterator itElm = (*it)->triangles.begin();
                  itElm != (*it)->triangles.end(); ++itElm)
                part[partIndex++] = mapOfPartitionsTag;
              for(std::vector<MQuadrangle*>::iterator itElm = (*it)->quadrangles.begin();
                  itElm != (*it)->quadrangles.end(); ++itElm)
                part[partIndex++] = mapOfPartitionsTag;
              mapOfPartitionsTag++;
            }
          }
          subGraph.partition(part);
    
          std::vector< std::set<MElement*> > subBoundaryElements =
            subGraph.getBoundaryElements(mapOfPartitionsTag);
    
          for(unsigned int i = 0; i < mapOfPartitionsTag; i++){
            for(std::set<MElement*>::iterator it = subBoundaryElements[i].begin();
                it != subBoundaryElements[i].end(); ++it){
              for(int j = 0; j < (*it)->getNumEdges(); j++){
                edgeToElement[(*it)->getEdge(j)].push_back
                (std::pair<MElement*, std::vector<unsigned int> >
                 (*it, mapOfPartitions[i]));
    
              }
            }
          }
        }
    
        int numEdgeEntity = model->getMaxElementaryNumber(1);
        for(hashmapedge::const_iterator it = edgeToElement.begin();
            it != edgeToElement.end(); ++it){
          MEdge e = it->first;
    
          std::vector<unsigned int> partitions;
          getPartitionInVector(partitions, it->second);
          if(partitions.size() < 2) continue;
    
          MElement* reference = getReferenceElement(it->second);
          if(!reference) continue;
    
          partitionEdge* pe = assignPartitionBoundary
            (model, e, reference, partitions, pedges, elementToEntity, numEdgeEntity);
          if(pe){
            std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
            for(unsigned int i = 0; i < it->second.size(); i++){
              boundaryEntityAndRefElement.insert
                (std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
                                                it->second[i].first));
            }
    
            assignBrep(model, boundaryEntityAndRefElement, pe);
          }
        }
        edgeToElement.clear();
    
        edges = model->getEdges();
        dividedNonConnectedEntities(model, 1, regions, faces, edges, vertices);
        elementToEntity.clear();
        fillElementToEntity(model, elementToEntity, 1);
      }
    
      if (meshDim >= 1){ // Create partition vertices
        Msg::Info(" - Creating partition vertices");
        if (meshDim == 1){
          for(unsigned int i = 0; i < model->getNumPartitions(); i++){
            for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
                it != boundaryElements[i].end(); ++it){
              for(int j = 0; j < (*it)->getNumPrimaryVertices(); j++){
                vertexToElement[(*it)->getVertex(j)].push_back
                (std::pair<MElement*, std::vector<unsigned int> >
                 (*it, std::vector<unsigned int>(1,i+1)));
              }
            }
          }
        }
        else{
          Graph subGraph(model);
          MakeGraph(model, subGraph, 1);
          subGraph.createDualGraph(false);
          std::vector<unsigned int> part(subGraph.ne());
          int partIndex = 0;
    
          std::map<unsigned int, std::vector<unsigned int> > mapOfPartitions;
          unsigned int mapOfPartitionsTag = 0;
          for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
            if((*it)->geomType() == GEntity::PartitionCurve){
              std::vector<unsigned int> partitions =
                static_cast<partitionEdge*>(*it)->getPartitions();
              mapOfPartitions.insert(std::pair<unsigned int, std::vector<unsigned int> >
                                     (mapOfPartitionsTag, partitions));
              // Must absolutely be in the same order as in the MakeGraph function
              for(std::vector<MLine*>::iterator itElm = (*it)->lines.begin();
                  itElm != (*it)->lines.end(); ++itElm)
                part[partIndex++] = mapOfPartitionsTag;
              mapOfPartitionsTag++;
            }
          }
          subGraph.partition(part);
    
          std::vector< std::set<MElement*> > subBoundaryElements =
            subGraph.getBoundaryElements(mapOfPartitionsTag);
    
          for(unsigned int i = 0; i < mapOfPartitionsTag; i++){
            for(std::set<MElement*>::iterator it = subBoundaryElements[i].begin();
                it != subBoundaryElements[i].end(); ++it){
              for(int j = 0; j < (*it)->getNumPrimaryVertices(); j++){
                vertexToElement[(*it)->getVertex(j)].push_back
                (std::pair<MElement*, std::vector<unsigned int> >
                 (*it, mapOfPartitions[i]));
              }
            }
          }
        }
        int numVertexEntity = model->getMaxElementaryNumber(1);
        for(hashmapvertex::const_iterator it = vertexToElement.begin();
            it != vertexToElement.end(); ++it){
          MVertex *v = it->first;
    
          std::vector<unsigned int> partitions;
          getPartitionInVector(partitions, it->second);
          if(partitions.size() < 2) continue;
    
          MElement* reference = getReferenceElement(it->second);
          if(!reference) continue;
    
          partitionVertex* pv = assignPartitionBoundary
            (model, v, reference, partitions, pvertices, elementToEntity, numVertexEntity);
          if(pv){
            std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
            for(unsigned int i = 0; i < it->second.size(); i++)
              boundaryEntityAndRefElement.insert
                (std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
                                                it->second[i].first));
    
            assignBrep(model, boundaryEntityAndRefElement, pv);
          }
        }
        vertexToElement.clear();
    
        vertices = model->getVertices();
        dividedNonConnectedEntities(model, 0, regions, faces, edges, vertices);
      }
    }
    
    static void addPhysical(GModel *const model, int level,
                            GEntity *childEntity, hashmap<std::string, int> &nameToNumber,
                            std::vector<GModel::piter> &iterators, int &numPhysical)
    {
      if(!CTX::instance()->mesh.partitionCreatePhysicals) return;
    
      unsigned int numPartitions = 0;
      if(childEntity->dim() == 3){
        numPartitions = static_cast<partitionRegion*>(childEntity)->numPartitions();
      }
      else if(childEntity->dim() == 2){
        numPartitions = static_cast<partitionFace*>(childEntity)->numPartitions();
      }
      else if(childEntity->dim() == 1){
        numPartitions = static_cast<partitionEdge*>(childEntity)->numPartitions();
      }
      else if(childEntity->dim() == 0){
        numPartitions = static_cast<partitionVertex*>(childEntity)->numPartitions();
      }
    
      std::string name;
      if(level == 0){
        name = "_omega{";
      }
      else{
        if(numPartitions== 1){
          name = "_omicron{";
        }
        else{
          if(level == 1){
            name = "_sigma{";
          }
          else if(level == 2){
            name = "_tau{";
          }
          else if(level == 3){
            name = "_upsilon{";
          }
        }
      }
    
      for(unsigned int i = 0; i < numPartitions; i++){
        if(i > 0) name += ",";
        unsigned int partition = 0;
        if(childEntity->dim() == 3){
          partition = static_cast<partitionRegion*>(childEntity)->getPartition(i);
        }
        else if(childEntity->dim() == 2){
          partition = static_cast<partitionFace*>(childEntity)->getPartition(i);
        }
        else if(childEntity->dim() == 1){
          partition = static_cast<partitionEdge*>(childEntity)->getPartition(i);
        }
        else if(childEntity->dim() == 0){
          partition = static_cast<partitionVertex*>(childEntity)->getPartition(i);
        }
    #if __cplusplus >= 201103L
        name += std::to_string(partition);
    #else
        char intToChar[20];
        sprintf(intToChar, "%d", partition);
        name += intToChar;
    #endif
      }
      name += "}";
    
      int number = 0;
      hashmap<std::string, int>::iterator it = nameToNumber.find(name);
      if(it == nameToNumber.end()){
        number = ++numPhysical;
        iterators[childEntity->dim()] = model->setPhysicalName
          (iterators[childEntity->dim()], name, childEntity->dim(), number);
        nameToNumber.insert(std::pair<std::string, int>(name, number));
      }
      else{
        number = it->second;
      }
      childEntity->addPhysicalEntity(number);
    }
    
    // AssignPhysicalName
    static void AssignPhysicalName(GModel *model)
    {
      std::vector<GModel::piter> iterators;
      model->getInnerPhysicalNamesIterators(iterators);
      int numPhysical = model->getMaxPhysicalNumber(-1);
      hashmap<GEntity*, int> levels;
      hashmap<std::string, int> nameToNumber;
      // Loop over regions
      for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
        if((*it)->geomType() == GEntity::PartitionVolume){
          addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
          levels.insert(std::pair<GEntity*, int>(*it, 0));
    
          std::vector<GFace*> listFace = (*it)->faces();
          for(std::vector<GFace*>::iterator itF = listFace.begin(); itF != listFace.end(); ++itF){
            if(levels.find(*itF) == levels.end()){
              const int level = levels[*it]+1;
              addPhysical(model, level, *itF, nameToNumber, iterators, numPhysical);
              levels.insert(std::pair<GEntity*, int>(*itF, level));
            }
          }
        }
      }
    
      // Loop over faces
      for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
        if((*it)->geomType() == GEntity::PartitionSurface){
          if(levels.find(*it) == levels.end()){
            addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
            levels.insert(std::pair<GEntity*, int>(*it, 0));
          }
    
          std::vector<GEdge*> const& listEdge = (*it)->edges();
          for(std::vector<GEdge*>::const_iterator itE = listEdge.begin(); itE != listEdge.end(); ++itE){
            if(levels.find(*itE) == levels.end()){
              const int level = levels[*it]+1;
              addPhysical(model, level, *itE, nameToNumber, iterators, numPhysical);
              levels.insert(std::pair<GEntity*, int>(*itE, level));
            }
          }
        }
      }
    
      // Loop over edges
      for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
        if((*it)->geomType() == GEntity::PartitionCurve){
          if(levels.find(*it) == levels.end()){
            addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
            levels.insert(std::pair<GEntity*, int>(*it, 0));
          }
    
          GVertex *v0 = (*it)->getBeginVertex();
          if(v0 && levels.find(v0) == levels.end()){
            const int level = levels[*it]+1;
            addPhysical(model, level, v0, nameToNumber, iterators, numPhysical);
            levels.insert(std::pair<GEntity*, int>(v0, level));
          }
          GVertex *v1 = (*it)->getEndVertex();
          if(v1 && levels.find(v1) == levels.end()){
            const int level = levels[*it]+1;
            addPhysical(model, level, v1, nameToNumber, iterators, numPhysical);
            levels.insert(std::pair<GEntity*, int>(v1, level));
          }
        }
      }
    }
    
    // Partition a mesh into n parts. Returns: 0 = success, 1 = error
    int PartitionMesh(GModel *const model)
    {
      if(CTX::instance()->mesh.numPartitions <= 0) return 0;
    
      Msg::StatusBar(true, "Partitioning mesh...");
      double t1 = Cpu();
    
      Graph graph(model);
      if(MakeGraph(model, graph, -1)) return 1;
      graph.nparts(CTX::instance()->mesh.numPartitions);
      if(PartitionGraph(graph)) return 1;
    
      // Assign partitions to elements
      hashmap<MElement*, unsigned int> elmToPartition;
      for(unsigned int i = 0; i < graph.ne(); i++){
        if(graph.element(i)){
          if(graph.nparts() > 1){
            elmToPartition.insert(std::pair<MElement*, unsigned int>(graph.element(i),
                                                                     graph.partition(i)+1));
            // Should be removed
            graph.element(i)->setPartition(graph.partition(i)+1);
          }
          else{
            elmToPartition.insert(std::pair<MElement*, unsigned int>(graph.element(i), 1));
            // Should be removed
            graph.element(i)->setPartition(1);
          }
        }
      }
      model->setNumPartitions(graph.nparts());
    
      CreateNewEntities(model, elmToPartition);
      elmToPartition.clear();
    
      double t2 = Cpu();
      Msg::StatusBar(true, "Done partitioning mesh (%g s)", t2 - t1);
    
      if(CTX::instance()->mesh.partitionCreateTopology){
        Msg::StatusBar(true, "Creating partition topology...");
        std::vector< std::set<MElement*> > boundaryElements = graph.getBoundaryElements();
        CreatePartitionTopology(model, boundaryElements, graph);
        boundaryElements.clear();
        AssignPhysicalName(model);
    
        double t3 = Cpu();
        Msg::StatusBar(true, "Done creating partition topology (%g s)", t3 - t2);
      }
    
      AssignMeshVertices(model);
    
      if(CTX::instance()->mesh.partitionCreateGhostCells){
        graph.clearDualGraph();
        graph.createDualGraph(true);
        graph.assignGhostCells();
      }
    
      return 0;
    }
    
    template <class ITERATOR>
    static void assignToParent(std::set<MVertex*> &verts, partitionRegion *region,
                               ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        if(region->getParentEntity()->dim() == 3)
          region->getParentEntity()->addElement((*it)->getType(), *it);
        (*it)->setPartition(0);
    
        for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
          if(verts.find((*it)->getVertex(i)) == verts.end()){
            (*it)->getVertex(i)->setEntity(region->getParentEntity());
            region->getParentEntity()->addMeshVertex((*it)->getVertex(i));
            verts.insert((*it)->getVertex(i));
          }
        }
      }
    }
    
    template <class ITERATOR>
    static void assignToParent(std::set<MVertex*> &verts, partitionFace *face,
                               ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        if(face->getParentEntity()->dim() == 2)
          face->getParentEntity()->addElement((*it)->getType(), *it);
        (*it)->setPartition(0);
    
        for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
          if(verts.find((*it)->getVertex(i)) == verts.end()){
            (*it)->getVertex(i)->setEntity(face->getParentEntity());
            face->getParentEntity()->addMeshVertex((*it)->getVertex(i));
            verts.insert((*it)->getVertex(i));
          }
        }
      }
    }
    
    template <class ITERATOR>
    static void assignToParent(std::set<MVertex*> &verts, partitionEdge *edge,
                               ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it){
        if(edge->getParentEntity()->dim() == 1)
          edge->getParentEntity()->addElement((*it)->getType(), *it);
        (*it)->setPartition(0);
    
        for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
          if(verts.find((*it)->getVertex(i)) == verts.end()){
            (*it)->getVertex(i)->setEntity(edge->getParentEntity());
            edge->getParentEntity()->addMeshVertex((*it)->getVertex(i));
            verts.insert((*it)->getVertex(i));
          }
        }
      }
    }
    
    template <class ITERATOR>
    void assignToParent(std::set<MVertex*> &verts, partitionVertex *vertex,
                        ITERATOR it_beg, ITERATOR it_end)
    {
      for(ITERATOR it = it_beg; it != it_end; ++it)
      {
        if(vertex->getParentEntity()->dim() == 0)
          vertex->getParentEntity()->addElement((*it)->getType(), *it);
        (*it)->setPartition(0);
    
        for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
          if(verts.find((*it)->getVertex(i)) == verts.end()){
            (*it)->getVertex(i)->setEntity(vertex->getParentEntity());
            vertex->getParentEntity()->addMeshVertex((*it)->getVertex(i));
            verts.insert((*it)->getVertex(i));
          }
        }
      }
    }
    
    // Un-partition a mesh and return to the initial mesh geomerty. Returns: 0 =
    // success, 1 = error.
    int UnpartitionMesh(GModel *const model)
    {
      // make a copy so we can iterate safely (we will remove some entities)
      std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
      std::set<GFace*, GEntityLessThan> faces = model->getFaces();
      std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
      std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
    
      std::set<MVertex*> verts;
    
      // Loop over vertices
      for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it){
        GVertex *vertex = *it;
    
        if(vertex->geomType() == GEntity::PartitionVertex){
          partitionVertex* pvertex = static_cast<partitionVertex*>(vertex);
          if(pvertex->getParentEntity()){
            assignToParent(verts, pvertex, pvertex->points.begin(), pvertex->points.end());
          }
          else{
            for(unsigned int j = 0; j < pvertex->points.size(); j++)
              delete pvertex->points[j];
          }
          pvertex->points.clear();
          pvertex->mesh_vertices.clear();
    
          model->remove(pvertex);
          delete pvertex;
        }
      }
    
      // Loop over edges
      for(GModel::eiter it = edges.begin(); it != edges.end(); ++it){
        GEdge *edge = *it;
        if(edge->geomType() == GEntity::PartitionCurve){
          partitionEdge* pedge = static_cast<partitionEdge*>(edge);
          if(pedge->getParentEntity()){
            assignToParent(verts, pedge, pedge->lines.begin(), pedge->lines.end());
          }
          else{
            for(unsigned int j = 0; j < pedge->lines.size(); j++)
              delete pedge->lines[j];
          }
          pedge->lines.clear();
          pedge->mesh_vertices.clear();
          pedge->setBeginVertex(0);
          pedge->setEndVertex(0);
    
          model->remove(pedge);
          delete pedge;
        }
        else if(edge->geomType() == GEntity::GhostCurve){
          model->remove(edge);
          delete edge;
        }
      }
    
      // Loop over faces
      for(GModel::fiter it = faces.begin(); it != faces.end(); ++it)
      {
        GFace *face = *it;
    
        if(face->geomType() == GEntity::PartitionSurface){
          partitionFace* pface = static_cast<partitionFace*>(face);
          if(pface->getParentEntity()){
            assignToParent(verts, pface, pface->triangles.begin(), pface->triangles.end());
            assignToParent(verts, pface, pface->quadrangles.begin(), pface->quadrangles.end());
          }
          else{
            for(unsigned int j = 0; j < pface->triangles.size(); j++)
              delete pface->triangles[j];
            for(unsigned int j = 0; j < pface->quadrangles.size(); j++)
              delete pface->quadrangles[j];
          }
          pface->triangles.clear();
          pface->quadrangles.clear();
          pface->mesh_vertices.clear();
          pface->set(std::vector<GEdge*>());
          pface->setOrientations(std::vector<int>());
    
          model->remove(pface);
          delete pface;
        }
        else if(face->geomType() == GEntity::GhostSurface){
          model->remove(face);
          delete face;
        }
      }
    
      // Loop over regions
      for(GModel::riter it = regions.begin(); it != regions.end(); ++it){
        GRegion *region = *it;
    
        if(region->geomType() == GEntity::PartitionVolume){
          partitionRegion* pregion = static_cast<partitionRegion*>(region);
          if(pregion->getParentEntity()){
            assignToParent(verts, pregion, pregion->tetrahedra.begin(), pregion->tetrahedra.end());
            assignToParent(verts, pregion, pregion->hexahedra.begin(), pregion->hexahedra.end());
            assignToParent(verts, pregion, pregion->prisms.begin(), pregion->prisms.end());
            assignToParent(verts, pregion, pregion->pyramids.begin(), pregion->pyramids.end());
            assignToParent(verts, pregion, pregion->trihedra.begin(), pregion->trihedra.end());
          }
          else{
            for(unsigned int j = 0; j < pregion->tetrahedra.size(); j++)
              delete pregion->tetrahedra[j];
    
            for(unsigned int j = 0; j < pregion->hexahedra.size(); j++)
              delete pregion->hexahedra[j];
    
            for(unsigned int j = 0; j < pregion->prisms.size(); j++)
              delete pregion->prisms[j];
    
            for(unsigned int j = 0; j < pregion->pyramids.size(); j++)
              delete pregion->pyramids[j];
    
            for(unsigned int j = 0; j < pregion->trihedra.size(); j++)
              delete pregion->trihedra[j];
          }
          pregion->tetrahedra.clear();
          pregion->hexahedra.clear();
          pregion->prisms.clear();
          pregion->pyramids.clear();
          pregion->trihedra.clear();
          pregion->mesh_vertices.clear();
          pregion->set(std::vector<GFace*>());
          pregion->setOrientations(std::vector<int>());
    
          model->remove(pregion);
          delete pregion;
        }
        else if(region->geomType() == GEntity::GhostVolume){
          model->remove(region);
          delete region;
        }
      }
    
      model->setNumPartitions(0);
    
      std::map<std::pair<int, int>, std::string> physicalNames = model->getPhysicalNames();
      for(GModel::piter it = physicalNames.begin(); it != physicalNames.end(); ++it){
        std::string name = it->second;
    
        if(name[0] == '_'){
          model->removePhysicalGroup(it->first.first, it->first.second);
        }
      }
    
      return 0;
    }
    
    // Create the partition according to the element split given by elmToPartition
    // Returns: 0 = success, 1 = no elements found.
    int PartitionUsingThisSplit(GModel *const model, unsigned int npart,
                                hashmap<MElement*, unsigned int> &elmToPartition)
    {
      Graph graph(model);
      if(MakeGraph(model, graph, -1)) return 1;
      graph.createDualGraph(false);
      graph.nparts(npart);
    
      if(elmToPartition.size() != graph.ne()){
        Msg::Error("All elements are not partitioned.");
        return 1;
      }
    
      std::vector<unsigned int> part(graph.ne());
      for(unsigned int i = 0; i < graph.ne(); i++){
        if(graph.element(i)){
          part[i] = elmToPartition[graph.element(i)]-1;
        }
      }
    
      // Check and correct the topology
      // Check and correct the topology
      for(unsigned int i = 0; i < graph.ne(); i++){
        if(graph.element(i)->getDim() == (int)graph.dim()) continue;
    
        for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
          if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()+1){
            if(part[i] != part[graph.adjncy(j)]){
              part[i] = part[graph.adjncy(j)];
              elmToPartition[graph.element(i)] = part[i]+1;
              break;
            }
          }
        }
    
        for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
          if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()+2){
            if(part[i] != part[graph.adjncy(j)]){
              part[i] = part[graph.adjncy(j)];
              elmToPartition[graph.element(i)] = part[i]+1;
              break;
            }
          }
        }
    
        for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
          if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()+3){
            if(part[i] != part[graph.adjncy(j)]){
              part[i] = part[graph.adjncy(j)];
              elmToPartition[graph.element(i)] = part[i]+1;
              break;
            }
          }
        }
      }
      graph.partition(part);
    
      model->setNumPartitions(graph.nparts());
    
      CreateNewEntities(model, elmToPartition);
      elmToPartition.clear();
    
      if(CTX::instance()->mesh.partitionCreateTopology){
        Msg::StatusBar(true, "Creating partition topology...");
        std::vector< std::set<MElement*> > boundaryElements = graph.getBoundaryElements();
        CreatePartitionTopology(model, boundaryElements, graph);
        boundaryElements.clear();
        AssignPhysicalName(model);
        Msg::StatusBar(true, "Done creating partition topology");
      }
    
      AssignMeshVertices(model);
    
      if(CTX::instance()->mesh.partitionCreateGhostCells){
        graph.clearDualGraph();
        graph.createDualGraph(false);
        graph.assignGhostCells();
      }
    
      return 0;
    }
    
    // Import a mesh partitionned by a tag given to the element and create the
    // topology (omega, sigma, bndSigma, ...). Returns: 0 = success, 1 = no elements
    // found.
    int ConvertOldPartitioningToNewOne(GModel *const model)
    {
      Msg::StatusBar(true, "Converting old partitioning...");
    
      hashmap<MElement*, unsigned int> elmToPartition;
      std::set<unsigned int> partitions;
      std::vector<GEntity*> entities;
      model->getEntities(entities);
      for(unsigned int i = 0; i < entities.size(); i++){
        for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
          MElement *e = entities[i]->getMeshElement(i);
          elmToPartition.insert(std::pair<MElement*, unsigned int>(e, e->getPartition()));
          partitions.insert(e->getPartition());
        }
      }
    
      return PartitionUsingThisSplit(model, partitions.size(), elmToPartition);
    }
    
    #else
    
    int PartitionMesh(GModel *const model)
    {
      Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
      return 0;
    }
    
    int UnpartitionMesh(GModel *const model)
    {
      return 0;
    }
    
    int ConvertOldPartitioningToNewOne(GModel *const model)
    {
      return 0;
    }
    
    int PartitionUsingThisSplit(GModel *const model, unsigned int npart,
                                hashmap<MElement*, unsigned int> &elmToPartition)
    {
      return 0;
    }
    
    #endif