Skip to content
Snippets Groups Projects
Select Git revision
  • b61da056ae01ac1e408e0dfc765feb344409b0cb
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • 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
  • gmsh_2_8_6
26 results

Options.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    • Christophe Geuzaine's avatar
      b61da056
      · b61da056
      Christophe Geuzaine authored
      extended selection buffer to prepare for mesh element picking and the new
      mesh editing interface
      b61da056
      History
      Christophe Geuzaine authored
      extended selection buffer to prepare for mesh element picking and the new
      mesh editing interface
    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