added partitioning files

Partitioning/CMakeLists.txt

+cmake_minimum_required(VERSION 3.5)
+
+project(graphPartitioner)
+
+option(MPI "Compile the parallel version" OFF)
+
+file(GLOB FILES main.cpp
+		io.cpp
+		graph.cpp
+		topology.cpp)
+
+include_directories(include/gmsh)
+include_directories(include/metis)
+
+if(MPI)
+	message("Parallel version")
+
+	set(CMAKE_CXX_COMPILER "mpic++")
+	set(CMAKE_C_COMPILER "mpicc")
+
+	find_package (MPI)
+	if (MPI_FOUND)
+		include_directories(${MPI_INCLUDE_PATH})
+	endif(MPI_FOUND)
+
+	add_definitions(-DPARALLEL)
+
+	file(GLOB FILESMPI main_mpi.cpp
+			   topology_mpi.cpp)
+else(MPI)
+	message("Sequential version")
+endif(MPI)
+
+if(APPLE)
+	message("OS X build")
+	link_directories(lib/osx)
+	set(EXECUTABLE_OUTPUT_PATH ../bin/osx)
+else(APPLE)
+	message("NIC4 build")
+	link_directories(lib/nic4)
+	set(EXECUTABLE_OUTPUT_PATH ../bin/nic4)
+endif(APPLE)
+
+
+if(MPI)
+	add_executable(meshPartitioner_mpi ${FILES} ${FILESMPI})
+	target_link_libraries(meshPartitioner_mpi Gmsh lapack blas parmetis metis)
+else(MPI)
+	add_executable(meshPartitioner ${FILES})
+	target_link_libraries(meshPartitioner Gmsh lapack blas metis)
+endif(MPI)

Partitioning/graph.cpp

+#include <map>
+#ifdef PARALLEL
+    #include <mpi.h>
+#endif
+
+#include "MElementCut.h"
+
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
+#include "MTrihedron.h"
+
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+
+#include "MLine.h"
+
+#include "MPoint.h"
+
+#include "graph.h"
+
+void SEQ::GModelToGraph(GModel* gModel, int* eptr, int** eind, int *metisToGmshIndex)
+{
+    std::multimap<int, int> elementsToNodes;
+
+    //Loop over regions
+    for(GModel::riter it = gModel->firstRegion(); it != gModel->lastRegion(); ++it)
+    {
+        GRegion *r = *it;
+    
+        fillElementsToNodes(elementsToNodes, r->tetrahedra.begin(), r->tetrahedra.end());
+        fillElementsToNodes(elementsToNodes, r->hexahedra.begin(), r->hexahedra.end());
+        fillElementsToNodes(elementsToNodes, r->prisms.begin(), r->prisms.end());
+        fillElementsToNodes(elementsToNodes, r->pyramids.begin(), r->pyramids.end());
+        fillElementsToNodes(elementsToNodes, r->trihedra.begin(), r->trihedra.end());
+        fillElementsToNodes(elementsToNodes, r->polyhedra.begin(), r->polyhedra.end());
+    }
+
+    //Loop over faces
+    for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
+    {
+        GFace *f = *it;
+
+        fillElementsToNodes(elementsToNodes, f->triangles.begin(), f->triangles.end());
+        fillElementsToNodes(elementsToNodes, f->quadrangles.begin(), f->quadrangles.end());
+        fillElementsToNodes(elementsToNodes, f->polygons.begin(), f->polygons.end());
+    }
+
+    //Loop over edges
+    for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
+    {
+        GEdge *e = *it;
+
+        fillElementsToNodes(elementsToNodes, e->lines.begin(), e->lines.end());
+    }
+
+    //Loop over vertices
+    for(GModel::viter it = gModel->firstVertex(); it != gModel->lastVertex(); ++it)
+    {
+        GVertex *v = *it;
+
+        fillElementsToNodes(elementsToNodes, v->points.begin(), v->points.end());
+    }
+
+    //create mesh format for METIS
+    unsigned int position = 0;
+    eptr[0] = 0;
+    unsigned int i = 0;
+    int itFirstLast = 0;
+    
+    for(std::multimap<int, int>::const_iterator it = elementsToNodes.begin(); it != elementsToNodes.end(); ++it)
+    {
+        if(itFirstLast != it->first+1)
+        {
+            metisToGmshIndex[i] = it->first + 1;
+            position += elementsToNodes.count(it->first);
+            eptr[i+1] = position;
+            i++;
+            itFirstLast = it->first+1;
+        }
+    }
+
+    (*eind) = new int[position];
+    unsigned int positionInd = 0;
+    for(std::multimap<int, int>::const_iterator it = elementsToNodes.begin(); it != elementsToNodes.end(); ++it)
+    {
+        int tag = it->second;
+        (*eind)[positionInd] = tag;
+        positionInd++;
+    }
+}
+
+template <class ITERATOR>
+void SEQ::fillElementsToNodes(std::multimap<int, int> &elementsToNodes, ITERATOR it_beg, ITERATOR it_end)
+{
+    for(ITERATOR it = it_beg; it != it_end; ++it)
+    {
+        const int tag = (*it)->getNum()-1;
+        
+        for(unsigned int j = 0; j < (*it)->getNumVertices(); j++)
+        {
+            elementsToNodes.insert(std::pair<int, int>(tag,(*it)->getVertex(j)->getNum()-1));
+        }
+    }
+}
+
+// Parallel
+
+#ifdef PARALLEL
+int PAR::GModelToGraph(GModel* gModel, int* eptr, int** eind, int *metisToGmshIndex, int* elmdist)
+{
+    int nbproc, myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &nbproc);
+    
+    std::multimap<int, int> elementsToNodes;
+    
+    for(unsigned int i = elmdist[myrank]; i < elmdist[myrank+1]; i++)
+    {
+        MElement* elm = gModel->getMeshElementByTag(i+1);
+        const int tag = i;
+        
+        for(unsigned int j = 0; j < elm->getNumVertices(); j++)
+        {
+            elementsToNodes.insert(std::pair<int, int>(tag,elm->getVertex(j)->getNum()-1));
+        }
+    }
+    
+    //create mesh format for METIS
+    unsigned int position = 0;
+    eptr[0] = 0;
+    unsigned int i = 0;
+    int itFirstLast = 0;
+    
+    for(std::multimap<int, int>::const_iterator it = elementsToNodes.begin(); it != elementsToNodes.end(); ++it)
+    {
+        if(itFirstLast != it->first+1)
+        {
+            metisToGmshIndex[i] = it->first + 1;
+            position += elementsToNodes.count(it->first);
+            eptr[i+1] = position;
+            i++;
+            itFirstLast = it->first+1;
+        }
+    }
+    
+    (*eind) = new int[position];
+    unsigned int positionInd = 0;
+    for(std::multimap<int, int>::const_iterator it = elementsToNodes.begin(); it != elementsToNodes.end(); ++it)
+    {
+        int tag = it->second;
+        (*eind)[positionInd] = tag;
+        positionInd++;
+    }
+    
+    return position;
+}
+#endif
+

Partitioning/graph.h

+#ifndef GRAPH_INCLUDED
+#define GRAPH_INCLUDED
+
+#include "Gmsh.h"
+#include "GModel.h"
+
+namespace SEQ {
+    void GModelToGraph(GModel* gModel, int* eptr, int** eind, int *metisToGmshIndex);
+    template <class ITERATOR>
+    void fillElementsToNodes(std::multimap<int, int> &nodesToElements, ITERATOR it_beg, ITERATOR it_end);
+}
+
+#ifdef PARALLEL
+namespace PAR {
+    int GModelToGraph(GModel* gModel, int* eptr, int** eind, int *metisToGmshIndex, int* elmdist);
+}
+#endif
+
+#endif //GRAPH_INCLUDED

Partitioning/io.cpp

+#include <map>
+#include <unordered_map>
+#include <vector>
+#include <fstream>
+
+#include "GModel.h"
+#include "GEntity.h"
+
+#include "partitionFace.h"
+#include "partitionEdge.h"
+#include "partitionVertex.h"
+
+#include "io.h"
+
+void SEQ::writeModels(std::vector<GModel*> &models)
+{
+    for(unsigned int i = 0; i < models.size(); i++)
+    {
+        std::string str = "mesh_";
+        str += std::to_string(i);
+        str += ".msh";
+        
+        models[i]->writeMSH(str.c_str());//, 3.0);
+    }
+}
+
+void SEQ::writeProFile(GModel* m, const int npart)
+{
+    std::ofstream file("partition.pro", std::ofstream::trunc);
+    
+    //-----------Group-----------
+    file << "Group{" << std::endl;
+    
+    //Omega
+    std::unordered_map<int, std::vector<int> > listOfOmega;//map between tag of omega and the physical's numbers that corresponds
+    for(GModel::piter it = m->firstPhysicalName(); it != m->lastPhysicalName(); ++it)
+    {
+        std::string name = it->second;
+        
+        if(name[0] == '_' && name[1] == 'o')
+        {
+            std::vector<int> num = getNumFromString(name);
+            
+            std::vector<int> vec = listOfOmega[num[0]];
+            vec.push_back(it->first.second);
+            listOfOmega[num[0]] = vec;
+        }
+    }
+    //Omega_i
+    for(std::unordered_map<int, std::vector<int> >::iterator it = listOfOmega.begin(); it != listOfOmega.end(); ++it)
+    {
+        std::vector<int> vec = it->second;
+        file << "\tOmega_" << it->first << " = Region[{";
+        
+        for(unsigned int i = 0; i < vec.size(); i++)
+        {
+            if(i != 0)
+            {
+                file << ", ";
+            }
+            file << vec[i];
+        }
+        file << "}];" << std::endl;
+    }
+    file << std::endl;
+    
+    //Sigma
+    if(npart == 1)
+    {
+        file << "\tSigma_0_0 = Region[{}];" << std::endl;
+        file << "\tBndSigma_0_0 =  Region[{}];" << std::endl;
+        file << "\tSigma_0 =  Region[{}];" << std::endl;
+        file << "\tBndGammaInf_0_0 = Region[{}];" << std::endl;
+        file << "\tBndGammaD_0_0 = Region[{}];" << std::endl;
+        file << "\tBndGammaInf_0 = Region[{}];" << std::endl;
+        file << "\tBndGammaD_0 = Region[{}];" << std::endl;
+        file << "\tD_0() = {0};" << std::endl;
+    }
+    std::unordered_map<int, std::vector<int> > listOfSigma;//map between tag of sigma and the physical's numbers that corresponds
+    std::unordered_map<int, std::vector<int> > listOfBndSigma;//map between tag of sigma's boundary and the physical's numbers that corresponds
+    for(GModel::piter it = m->firstPhysicalName(); it != m->lastPhysicalName(); ++it)
+    {
+        std::string name = it->second;
+        
+        if(name[0] == '_' && name[1] == 's')
+        {
+            std::vector<int> num = getNumFromString(name);
+            
+            if(num.size() < 3)
+            {
+                for(unsigned int i = 0; i < num.size(); i++)
+                {
+                    std::vector<int> vec = listOfSigma[num[i]];
+                    vec.push_back(it->first.second);
+                    listOfSigma[num[i]] = vec;
+                }
+            }
+            else
+            {
+                for(unsigned int i = 0; i < num.size(); i++)
+                {
+                    std::vector<int> vec = listOfBndSigma[num[i]];
+                    vec.push_back(it->first.second);
+                    listOfBndSigma[num[i]] = vec;
+                }
+            }
+        }
+    }
+    file << std::endl;
+    //Sigma_i_j and BndSigma_i_j
+    std::unordered_map<int, std::vector<int> > listOfNeighbour;//map between tag of omega and tag of neighbours
+    for(std::unordered_map<int, std::vector<int> >::iterator it = listOfSigma.begin(); it != listOfSigma.end(); ++it)
+    {
+        for(std::unordered_map<int, std::vector<int> >::iterator it2 = it; it2 != listOfSigma.end(); ++it2)
+        {
+            if(it != it2)
+            {
+                std::vector<int> vec1 = it->second;
+                std::vector<int> vec2 = it2->second;
+                std::vector<int>* vecCommun =  new std::vector<int>;
+            
+                if(communPhysicals(vec1, vec2, vecCommun))
+                {
+                    listOfNeighbour[it->first].push_back(it2->first);
+                    listOfNeighbour[it2->first].push_back(it->first);
+                    
+                    file << "\tSigma_" << it->first << "_" << it2->first << " = Region[{";
+                    for(unsigned int i = 0; i < vecCommun->size(); i++)
+                    {
+                        if(i != 0)
+                        {
+                            file << ", ";
+                        }
+                        file << (*vecCommun)[i];
+                    }
+                    file << "}];" << std::endl;
+                    
+                    file << "\tSigma_" << it2->first << "_" << it->first << " = Region[{";
+                    for(unsigned int i = 0; i < vecCommun->size(); i++)
+                    {
+                        if(i != 0)
+                        {
+                            file << ", ";
+                        }
+                        file << (*vecCommun)[i];
+                    }
+                    file << "}];" << std::endl;
+                    
+                    if(listOfBndSigma.count(it->first) > 0)
+                    {
+                        std::vector<int> vec1 = listOfBndSigma[it->first];
+                        std::vector<int> vec2 = listOfBndSigma[it2->first];
+                        std::vector<int>* vecCommun =  new std::vector<int>;
+                        
+                        if(communPhysicals(vec1, vec2, vecCommun))
+                        {
+                            file << "\tBndSigma_" << it->first << "_" << it2->first << " = Region[{";
+                            for(unsigned int i = 0; i < vecCommun->size(); i++)
+                            {
+                                if(i != 0)
+                                {
+                                    file << ", ";
+                                }
+                                file << (*vecCommun)[i];
+                            }
+                            file << "}];" << std::endl;
+                            
+                            file << "\tBndSigma_" << it2->first << "_" << it->first << " = Region[{";
+                            for(unsigned int i = 0; i < vecCommun->size(); i++)
+                            {
+                                if(i != 0)
+                                {
+                                    file << ", ";
+                                }
+                                file << (*vecCommun)[i];
+                            }
+                            file << "}];" << std::endl;
+                        }
+                        else
+                        {
+                            file << "\tBndSigma_" << it->first << "_" << it2->first << " = Region[{}];" << std::endl;
+                            file << "\tBndSigma_" << it2->first << "_" << it->first << " = Region[{}];" << std::endl;
+                        }
+                    }
+                    else
+                    {
+                        file << "\tBndSigma_" << it->first << "_" << it2->first << " = Region[{}];" << std::endl;
+                        file << "\tBndSigma_" << it2->first << "_" << it->first << " = Region[{}];" << std::endl;
+                    }
+                    
+                    file << "\tBndGammaInf_" << it->first << "_" << it2->first << " = Region[{}];" << std::endl;
+                    file << "\tBndGammaInf_" << it2->first << "_" << it->first << " = Region[{}];" << std::endl;
+                    
+                    file << "\tBndGammaD_" << it->first << "_" << it2->first << " = Region[{}];" << std::endl;
+                    file << "\tBndGammaD_" << it2->first << "_" << it->first << " = Region[{}];" << std::endl;
+                    
+                    file << "\tBndGammaInf_" << it->first << " = Region[{}];" << std::endl;
+                    file << "\tBndGammaInf_" << it2->first << " = Region[{}];" << std::endl;
+                    
+                    file << "\tBndGammaD_" << it->first << " = Region[{}];" << std::endl;
+                    file << "\tBndGammaD_" << it2->first << " = Region[{}];" << std::endl;
+                }
+                delete vecCommun;
+            }
+        }
+    }
+    file << std::endl;
+    //Sigma_i
+    for(std::unordered_map<int, std::vector<int> >::iterator it = listOfSigma.begin(); it != listOfSigma.end(); ++it)
+    {
+        std::vector<int> vec = it->second;
+        file << "\tSigma_" << it->first << " = Region[{";
+        
+        for(unsigned int i = 0; i < vec.size(); i++)
+        {
+            if(i != 0)
+            {
+                file << ", ";
+            }
+            file << vec[i];
+        }
+        file << "}];" << std::endl;
+    }
+    file << std::endl;
+    //BndSigma_i
+    for(std::unordered_map<int, std::vector<int> >::iterator it = listOfBndSigma.begin(); it != listOfBndSigma.end(); ++it)
+    {
+        std::vector<int> vec = it->second;
+        file << "\tBndSigma_" << it->first << " = Region[{";
+        
+        for(unsigned int i = 0; i < vec.size(); i++)
+        {
+            if(i != 0)
+            {
+                file << ", ";
+            }
+            file << vec[i];
+        }
+        file << "}];" << std::endl;
+    }
+    file << std::endl << std::endl;
+    
+    //D
+    file << "\tD() = {";
+    for(unsigned int i = 0; i < listOfOmega.size(); i++)
+    {
+        if(i != 0)
+        {
+            file << ", ";
+        }
+        file << i;
+    }
+    file << "};" << std::endl;
+    file << "\tN_DOM = #D();" << std::endl;
+    
+    //D_i
+    for(std::unordered_map<int, std::vector<int> >::iterator it = listOfNeighbour.begin(); it != listOfNeighbour.end(); ++it)
+    {
+        file << "\tD_" << it->first << " = {";
+        for(unsigned int i = 0; i < it->second.size(); i++)
+        {
+            if(i != 0)
+            {
+                file << ", ";
+            }
+            file << it->second[i];
+        }
+        file << "};" << std::endl;
+    }
+    
+    
+    file << "}" << std::endl << std::endl << std::endl;
+    //**********************************************
+    
+    
+    //-----------Function-----------
+    file << "Function {" << std::endl;
+    file << "\tmyD = {} ; // the domains that I'm in charge of" << std::endl;
+    for(unsigned int i = 0; i < listOfOmega.size(); i++)
+    {
+        file << "\tmyD_" << i << " = {};" << std::endl;
+    }
+    file << "\tListOfFields = {};" << std::endl;
+    file << "\tListOfConnectedFields = {};" << std::endl << std::endl;
+    
+    
+    file << "\tFor idom In {0:N_DOM-1}" << std::endl;
+    
+    file << "\t\tIf (idom % MPI_Size == MPI_Rank)" << std::endl;
+    
+    file << "\t\t\tmyD() += D(idom);" << std::endl;
+    file << "\t\t\tmyD~{idom} += D~{idom}();" << std::endl;
+    
+    file << "\t\tEndIf" << std::endl;
+    file << "\tEndFor" << std::endl;
+    
+    file << "\tFor ii In {0:#myD()-1}" << std::endl;
+    
+    file << "\t\ti = myD(ii);" << std::endl;
+    
+    file << "\t\tIf(#myD~{i}() == 2)" << std::endl;
+    file << "\t\t\tPrintf(\"We can do sweeping!\");" << std::endl;
+    file << "\t\tEndIf" << std::endl;
+
+    file << "\t\tFor jj In {0:#myD~{i}()-1}" << std::endl;
+    file << "\t\t\tj = myD~{i}(jj);" << std::endl << std::endl;
+    
+    file << "\t\t\ttag_g~{i}~{j} = i * 1000 + j;" << std::endl;
+    file << "\t\t\ttag_g~{j}~{i} = j * 1000 + i;" << std::endl << std::endl;
+    
+    file << "\t\t\tListOfFields() += tag_g~{i}~{j};" << std::endl;
+    file << "\t\t\tListOfConnectedFields() += 1;" << std::endl;
+    file << "\t\t\tListOfConnectedFields() += tag_g~{j}~{i};" << std::endl;
+    
+    file << "\t\t\tIf(ANALYSIS == 0)" << std::endl;
+    file << "\t\t\t\tg_in~{i}~{j}[ Sigma~{i}~{j} ] = ComplexScalarField[XYZ[]]{ tag_g~{j}~{i} };" << std::endl;
+    file << "\t\t\tEndIf" << std::endl;
+    file << "\t\t\tIf(ANALYSIS == 1)" << std::endl;
+    file << "\t\t\t\tg_in~{i}~{j}[ Sigma~{i}~{j} ] = ComplexVectorField[XYZ[]]{ tag_g~{j}~{i} };" << std::endl;
+    file << "\t\t\tEndIf" << std::endl;
+    
+    file << "\t\tEndFor" << std::endl;
+
+    file << "\tEndFor" << std::endl;
+    
+    file << "}" << std::endl;
+    
+    file.close();
+}
+
+std::vector<int> SEQ::getNumFromString(std::string name)
+{
+    std::vector<int> num;
+    std::string currentNum;
+    
+    for(unsigned int i = 0; i < name.size(); i++)
+    {
+        if(name[i] == '0' || name[i] == '1' || name[i] == '2'|| name[i] == '3'|| name[i] == '4'|| name[i] == '5'|| name[i] == '6'|| name[i] == '7'|| name[i] == '8'|| name[i] == '9')
+        {
+            currentNum += name[i];
+        }
+        
+        if(name[i] == ',' || name[i] == '}')
+        {
+            num.push_back(stoi(currentNum));
+            currentNum.clear();
+        }
+    }
+    
+    return num;
+}
+
+bool SEQ::communPhysicals(const std::vector<int> vec1, const std::vector<int> vec2, std::vector<int>* vecCommun)
+{
+    for(unsigned int i = 0; i < vec1.size(); i++)
+    {
+        for(unsigned int j = 0; j < vec2.size(); j++)
+        {
+            if(vec1[i] == vec2[j])
+            {
+                vecCommun->push_back(vec1[i]);
+            }
+        }
+    }
+    
+    if(vecCommun->size() > 0)
+    {
+        return true;
+    }
+    return false;
+}
+
+
+

Partitioning/io.h

+#ifndef IO_INCLUDED
+#define IO_INCLUDED
+
+namespace SEQ {
+    void writeModels(std::vector<GModel*> &models);
+    void writeProFile(GModel* m, const int npart);
+    std::vector<int> getNumFromString(std::string name);
+    bool communPhysicals(const std::vector<int> vec1, const std::vector<int> vec2, std::vector<int>* vecCommun);
+}
+
+#endif /* IO_INCLUDED */

Partitioning/main.cpp

+#include <iostream>
+#include <cstdlib>
+#include <string>
+#include <ctime>
+#include <fstream>
+
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+
+#ifdef PARALLEL
+    #include <mpi.h>
+    #include "main_mpi.h"
+#else
+    #include "metis.h"
+#endif
+
+#include "graph.h"
+#include "topology.h"
+#include "io.h"
+
+using namespace SEQ;
+
+int main(int argc, char **argv)
+{
+#ifdef PARALLEL
+    int nbproc, myrank;
+    MPI_Init(&argc,&argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &nbproc);
+    
+    clock_t t1, t2;
+    if(myrank == 0)
+    {
+        std::cout << "#############################################################################" << std::endl;
+        std::cout << "#                                                                           #" << std::endl;
+        std::cout << "#                         #   #  #####   ####  #   #                        #" << std::endl;
+        std::cout << "#                         ## ##  #      #      #   #                        #" << std::endl;
+        std::cout << "#                         # # #  ####   #####  #####                        #" << std::endl;
+        std::cout << "#                         #   #  #          #  #   #                        #" << std::endl;
+        std::cout << "#                         #   #  #####  ####   #   #                        #" << std::endl;
+        std::cout << "#                                                                           #" << std::endl;
+        std::cout << "#                                                                           #" << std::endl;
+        std::cout << "#  ####    ###   #####  #####  ###  #####  ###   ###   #   #  #####  #####  #" << std::endl;
+        std::cout << "#  #   #  #   #  #   #    #     #     #     #   #   #  ##  #  #      #   #  #" << std::endl;
+        std::cout << "#  ####   #####  #####    #     #     #     #   #   #  # # #  ####   #####  #" << std::endl;
+        std::cout << "#  #      #   #  # ##     #     #     #     #   #   #  # # #  #      # ##   #" << std::endl;
+        std::cout << "#  #      #   #  #   #    #    ###    #    ###   ###   #  ##  #####  #   #  #" << std::endl;
+        std::cout << "#                                                                           #" << std::endl;
+        std::cout << "####################################################### Parallel version ####" << std::endl << std::endl;
+        
+        if(argc < 2)
+        {
+            std::cout << "Arguments missing! Syntaxe :" << std::endl;
+            std::cout << argv[0] << " \"mesh.msh\" \"nbrPartitions\"" << std::endl;
+            MPI_Finalize();
+            return 0;
+        }
+        
+        if(atoi(argv[2]) < 1)
+        {
+            std::cout << "The number of partition must be greater than zero!" << std::endl;
+            MPI_Finalize();
+            return 0;
+        }
+        t1 = clock();
+    }
+    
+    MPI_Barrier(MPI_COMM_WORLD);
+    GmshInitialize(argc, argv);
+    main_mpi(nbproc, myrank, argv[1], atoi(argv[2]));
+    MPI_Barrier(MPI_COMM_WORLD);
+    GmshFinalize();
+    
+    if(myrank == 0)
+    {
+        t2 = clock();
+        float temps = (float)(t2-t1)/CLOCKS_PER_SEC;
+        
+        std::cout << "-> Partition done in " << temps << "seconds" << std::endl;
+        
+        std::ofstream time("time.txt", std::ofstream::app);
+        time << temps << std::endl;
+        time.close();
+    }
+    MPI_Finalize();
+#else
+    std::cout << "#############################################################################" << std::endl;
+    std::cout << "#                                                                           #" << std::endl;
+    std::cout << "#                         #   #  #####   ####  #   #                        #" << std::endl;
+    std::cout << "#                         ## ##  #      #      #   #                        #" << std::endl;
+    std::cout << "#                         # # #  ####   #####  #####                        #" << std::endl;
+    std::cout << "#                         #   #  #          #  #   #                        #" << std::endl;
+    std::cout << "#                         #   #  #####  ####   #   #                        #" << std::endl;
+    std::cout << "#                                                                           #" << std::endl;
+    std::cout << "#                                                                           #" << std::endl;
+    std::cout << "#  ####    ###   #####  #####  ###  #####  ###   ###   #   #  #####  #####  #" << std::endl;
+    std::cout << "#  #   #  #   #  #   #    #     #     #     #   #   #  ##  #  #      #   #  #" << std::endl;
+    std::cout << "#  ####   #####  #####    #     #     #     #   #   #  # # #  ####   #####  #" << std::endl;
+    std::cout << "#  #      #   #  # ##     #     #     #     #   #   #  # # #  #      # ##   #" << std::endl;
+    std::cout << "#  #      #   #  #   #    #    ###    #    ###   ###   #  ##  #####  #   #  #" << std::endl;
+    std::cout << "#                                                                           #" << std::endl;
+    std::cout << "#############################################################################" << std::endl << std::endl;
+    
+    if(argc < 2)
+    {
+        std::cout << "Arguments missing! Syntaxe :" << std::endl;
+        std::cout << argv[0] << " \"mesh.msh\" \"nbrPartitions\"" << std::endl;
+        return 0;
+    }
+    
+    if(atoi(argv[2]) < 1)
+    {
+        std::cout << "The number of partition must be greater than zero!" << std::endl;
+        return 0;
+    }
+    
+    float temps;
+    clock_t t1, t2;
+    t1 = clock();
+    
+    GmshInitialize(argc, argv);
+    GModel *m = new GModel();
+    std::cout << "Reading msh file... " << std::flush;
+    m->readMSH(argv[1]);
+    std::cout << "Done!" << std::endl;
+
+    const int numElements = m->getNumMeshElements();
+    const int numVertices = m->getNumMeshVertices();
+    
+    idx_t nparts = atoi(argv[2]);
+    
+    if(nparts > 1)
+    {
+        idx_t *epart = new idx_t[numElements];
+        idx_t *npart = new idx_t[numVertices];
+        int* eptr = new int[numElements+1];
+        int* eind = NULL;
+        int *metisToGmshIndex = new int[numElements];
+        
+        std::cout << "Creating Metis structure... " << std::flush;
+        GModelToGraph(m, eptr, &eind, metisToGmshIndex);
+        std::cout << "Done!" << std::endl;
+        
+        std::cout << "Mesh partitioning... " << std::flush;
+        idx_t objval;
+        idx_t ncommon = m->getDim();
+        idx_t options[METIS_NOPTIONS];
+        METIS_SetDefaultOptions(options);
+    
+        const int error = METIS_PartMeshDual((idx_t*)&numElements, (idx_t*)&numVertices, (idx_t*)eptr, (idx_t*)eind, NULL, NULL, &ncommon, &nparts, NULL, options, &objval, epart, npart);
+    
+        switch(error)
+        {
+            case METIS_OK:
+                std::cout << "Done!" << std::endl;
+                break;
+            case METIS_ERROR_INPUT:
+                std::cout << "Metis error (input)!" << std::endl;
+                return 0;
+                break;
+            case METIS_ERROR_MEMORY:
+                std::cout << "Metis error (memory)!" << std::endl;
+                return 0;
+                break;
+            case METIS_ERROR:
+                std::cout << "Metis error!" << std::endl;
+                return 0;
+                break;
+            default:
+                std::cout << "Error!" << std::endl;
+                return 0;
+                break;
+        }
+
+        for(unsigned int i = 0; i < numElements; i++)
+        {
+            m->getMeshElementByTag(metisToGmshIndex[i])->setPartition(epart[i]);
+        }
+        
+        delete[] eptr;
+        delete[] eind;
+        delete[] epart;
+        delete[] npart;
+        delete[] metisToGmshIndex;
+    }
+    
+    std::cout << "Creating new GModel... " << std::flush;
+    GModel* global = new GModel();
+    std::vector<GModel*> models = createNewModels(m, global, nparts);
+    std::cout << "Done!" << std::endl;
+    
+    std::cout << "Assign mesh vertices to models... " << std::flush;
+    for (unsigned int i = 0; i < nparts; i++)
+    {
+        assignMeshVerticesToModel(models[i]);
+    }
+    assignMeshVerticesToModel(global);
+    std::cout << "Done!" << std::endl;
+    
+    std::cout << "Creating new elements... " << std::endl;
+    createPartitionBoundaries(global, models);
+    std::cout << "Done!" << std::endl;
+    
+    std::cout << "Writing partition meshes... " << std::flush;
+    writeModels(models);
+    std::cout << "Done!" << std::endl;
+    
+    std::cout << "Writing global mesh... " << std::flush;
+    global->writeMSH("global.msh", 3.0);
+    std::cout << "Done!" << std::endl;
+    
+    std::cout << "Writing .pro file... " << std::flush;
+    writeProFile(global, nparts);
+    std::cout << "Done!" << std::endl;
+    
+    freeModels(models, global);
+
+    delete m;
+    GmshFinalize();
+    
+    t2 = clock();
+    temps = (float)(t2-t1)/CLOCKS_PER_SEC;
+    
+    std::cout << "-> Partition done in " << temps << "seconds" << std::endl;
+    
+    std::ofstream time("time.txt", std::ofstream::app);
+    time << temps << std::endl;
+    time.close();
+#endif
+    
+    return 0;
+}

Partitioning/main_mpi.cpp

+#include <iostream>
+#include <mpi.h>
+
+#include "parmetis.h"
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+
+#include "main_mpi.h"
+#include "graph.h"
+#include "topology.h"
+#include "topology_mpi.h"
+#include "io.h"
+
+using namespace PAR;
+
+int main_mpi(int nbproc, int myrank, std::string path, int nparts)
+{
+    GModel *m = new GModel();
+    if(myrank == 0) std::cout << "Reading msh file... " << std::flush;
+    m->readMSH(path.c_str());
+    if(myrank == 0) std::cout << "Done!" << std::endl;
+    
+    
+    const int numElements = m->getNumMeshElements();
+    const int numVertices = m->getNumMeshVertices();
+    
+    if(nparts > 1)
+    {
+        int *elmdist = new int[nbproc+1];
+        elmdist[0] = 0;
+        for(unsigned int i = 1; i < nbproc; i++)
+        {
+            elmdist[i] = i*numElements/nbproc;
+        }
+        elmdist[nbproc] = numElements;
+        
+        /*for(unsigned int i = 0; i < nbproc+1; i++)
+        {
+            if(myrank == 0) std::cout << elmdist[i] << std::endl;
+        }*/
+        
+        
+        int* eptr = new int[(elmdist[myrank+1]-elmdist[myrank])+1];
+        int* eind = NULL;
+        int *metisToGmshIndex = new int[elmdist[myrank+1]-elmdist[myrank]];
+        
+        if(myrank == 0) std::cout << "Creating Metis structure... " << std::flush;
+        const int pos = GModelToGraph(m, eptr, &eind, metisToGmshIndex, elmdist);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+        
+        
+        if(myrank == 0) std::cout << "Mesh partitioning... " << std::flush;
+        idx_t numflag = 0;
+        idx_t ncommon = m->getDim();
+        idx_t *xadj = NULL;
+        idx_t *adjncy = NULL;
+        MPI_Comm comm = MPI_COMM_WORLD;
+        
+        int error = ParMETIS_V3_Mesh2Dual((idx_t*)elmdist, (idx_t*)eptr, (idx_t*)eind, &numflag, &ncommon, &xadj, &adjncy, &comm);
+        switch(error)
+        {
+            case METIS_OK:
+                break;
+            case METIS_ERROR:
+                std::cout << "[rank " << myrank << "] Metis error!" << std::endl;
+                MPI_Finalize();
+                return 0;
+                break;
+        }
+        
+        idx_t objval;
+        idx_t wgtflag = 0;
+        idx_t ncon = 1;
+        real_t *tpwgts = new real_t[nparts];
+        for(unsigned int i = 0; i < nparts; i++)
+        {
+            tpwgts[i] = 1./nparts;
+        }
+        real_t ubvec[1];
+        ubvec[0] = 1.05;
+        idx_t options[3];
+        options[0] = 0;
+        idx_t* part = new idx_t[elmdist[myrank+1]-elmdist[myrank]];
+        error = ParMETIS_V3_PartKway((idx_t*)elmdist, xadj, adjncy, NULL, NULL, &wgtflag, &numflag, &ncon, (idx_t*)&nparts, tpwgts, ubvec, options, &objval, part, &comm);
+        
+        /*for(unsigned int i = 0; i < elmdist[myrank+1]-elmdist[myrank]; i++)
+        {
+            std::cout << "[rank " << myrank << "] " << part[i] << std::endl;
+        }*/
+        
+        switch(error)
+        {
+            case METIS_OK:
+                if(myrank == 0) std::cout << "Done!" << std::endl;
+                break;
+            case METIS_ERROR:
+                std::cout << "[rank " << myrank << "] Metis error!" << std::endl;
+                MPI_Finalize();
+                return 0;
+                break;
+        }
+        
+        MPI_Barrier(MPI_COMM_WORLD);
+        for(unsigned int i = 0; i < nbproc; i++)
+        {
+            if(i == myrank)
+            {
+                for(unsigned int j = 0; j < elmdist[myrank+1]-elmdist[myrank]; j++)
+                {
+                    m->getMeshElementByTag(metisToGmshIndex[j])->setPartition(part[j]);
+                    
+                    for(unsigned int k = 0; k < nbproc; k++)
+                    {
+                        if(k != myrank)
+                        {
+                            //element tag
+                            MPI_Send(&metisToGmshIndex[j], 1, MPI_INT, k, 1, MPI_COMM_WORLD);
+                            //partition
+                            MPI_Send(&part[j], 1, MPI_INT, k, 2, MPI_COMM_WORLD);
+                        }
+                    }
+                }
+            }
+            else
+            {
+                for(unsigned int j = 0; j < elmdist[i+1]-elmdist[i]; j++)
+                {
+                    int tag = 0, numPart = 0;
+                    MPI_Status status;
+                    //element tag
+                    MPI_Recv(&tag, 1, MPI_INT, i, 1, MPI_COMM_WORLD, &status);
+                    //partition
+                    MPI_Recv(&numPart, 1, MPI_INT, i, 2, MPI_COMM_WORLD, &status);
+                    
+                    m->getMeshElementByTag(tag)->setPartition(numPart);
+                }
+            }
+        }
+        
+        delete[] eptr;
+        delete[] eind;
+        delete[] metisToGmshIndex;
+        delete[] xadj;
+        delete[] adjncy;
+        delete[] tpwgts;
+        delete[] part;
+        delete[] elmdist;
+    }
+    
+    std::vector<int> myPart;
+    for(unsigned int i = 0; i < nparts; i++)
+    {
+        if((i+1)%nbproc == myrank)
+        {
+            myPart.push_back(i);
+        }
+    }
+    
+    if(myPart.size() != 0)
+    {
+        if(myrank == 0) std::cout << "Creating new GModel... " << std::flush;
+        GModel* global = new GModel();
+        std::vector<GModel*> models = createNewModels(m, global, nparts, myPart);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+            
+        if(myrank == 0) std::cout << "Assign mesh vertices to models... " << std::flush;
+        for (unsigned int i = 0; i < myPart.size(); i++)
+        {
+            SEQ::assignMeshVerticesToModel(models[i]);
+        }
+        if(myrank == 0) SEQ::assignMeshVerticesToModel(global);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+            
+        if(myrank == 0) std::cout << "Creating new elements... " << std::endl;
+        SEQ::createPartitionBoundaries(global, models);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+        
+        if(myrank == 0) std::cout << "Writing partition meshes... " << std::flush;
+        SEQ::writeModels(models);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+        
+        if(myrank == 0) std::cout << "Writing global mesh... " << std::flush;
+        if(myrank == 0) global->writeMSH("global.msh", 3.0);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+        
+        if(myrank == 0) std::cout << "Writing .pro file... " << std::flush;
+        if(myrank == 0) SEQ::writeProFile(global, nparts);
+        if(myrank == 0) std::cout << "Done!" << std::endl;
+        
+        SEQ::freeModels(models, global);
+    }
+    else
+    {
+        std::cout << "[rank " << myrank << "] Unused rank because nparts is smaller than nbproc!" << std::endl;
+    }
+    
+    delete m;
+    return 0;
+}

Partitioning/main_mpi.h

+#ifndef MAIN_MPI_INCLUDED
+#define MAIN_MPI_INCLUDED
+
+#include <string>
+
+int main_mpi(int nbproc, int myrank, std::string path, int nparts);
+
+#endif //MAIN_MPI_INCLUDED

Partitioning/topology.h

+#ifndef TOPOLOGY_INCLUDED
+#define TOPOLOGY_INCLUDED
+
+#include <unordered_map>
+
+#include "GModel.h"
+
+#include "partitionFace.h"
+#include "partitionEdge.h"
+#include "partitionVertex.h"
+
+#include "MFaceHash.h"
+#include "MEdgeHash.h"
+
+namespace SEQ {
+    int createPartitionBoundaries(GModel *model, std::vector<GModel*> models);
+    
+    template <class ITERATOR>
+    void fillit_(std::unordered_map<MFace, std::vector<MElement*> , Hash_Face, Equal_Face> &faceToElement, ITERATOR it_beg, ITERATOR it_end);
+    template <class ITERATOR>
+    void fillit_(std::unordered_map<MEdge, std::vector<MElement*> , Hash_Edge, Equal_Edge> &edgeToElement, ITERATOR it_beg, ITERATOR it_end);
+    template <class ITERATOR>
+    void fillit_(std::unordered_map<MVertex*, std::vector<MElement*> > &vertexToElement, ITERATOR it_beg, ITERATOR it_end);
+    
+    void assignPartitionBoundary(GModel *model, MFace &me, std::set<partitionFace*, Less_partitionFace> &pfaces, std::vector<MElement*> &v, std::vector<GModel*> models, std::vector< std::set<partitionFace*, Less_partitionFace> > &pfacesOfModels);
+    void assignPartitionBoundary(GModel *model, MEdge &me, std::set<partitionEdge*, Less_partitionEdge> &pedges, std::vector<MElement*> &v, std::set<partitionFace*, Less_partitionFace> &pfaces, std::vector<GModel*> models, std::vector< std::set<partitionEdge*, Less_partitionEdge> > &pedgesOfModels);
+    void assignPartitionBoundary(GModel *model, MVertex *ve, std::set<partitionVertex*, Less_partitionVertex> &pvertices, std::vector<MElement*> &v, std::set<partitionEdge*, Less_partitionEdge> &pedges, std::set<partitionFace*, Less_partitionFace> &pfaces, std::vector<GModel*> models, std::vector< std::set<partitionVertex*, Less_partitionVertex> > &pverticesOfModels);
+    
+    std::vector<GModel*> createNewModels(GModel *gModel, GModel *global, int nparts);
+    
+    void assignMeshVerticesToModel(GModel *gModel);
+    template <class ITERATOR>
+    void fillVertexToEntity(std::unordered_map<MVertex*, GEntity*> &vertexToEntity, GEntity* entity, ITERATOR it_beg, ITERATOR it_end);
+        
+    void freeModels(std::vector<GModel*> &models, GModel *global);
+    
+    void addPhysical(GModel *newModel, GEntity *newEntity, GModel *oldModel, GEntity *oldEntity, GModel *globalModel, GEntity *globalEntity, int partition, int maxDim);
+}
+
+
+#endif //TOPOLOGY_INCLUDED

Partitioning/topology_mpi.cpp

+#include <iostream>
+#include <set>
+#include <map>
+#include <unordered_map>
+#include <string>
+#include <algorithm>
+#include <mpi.h>
+
+#include "MElement.h"
+#include "MElementCut.h"
+
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
+#include "MTrihedron.h"
+
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+
+#include "MPoint.h"
+
+#include "MElementCut.h"
+
+#include "discreteRegion.h"
+#include "discreteFace.h"
+#include "discreteEdge.h"
+#include "discreteVertex.h"
+
+#include "partitionFace.h"
+#include "partitionEdge.h"
+#include "partitionVertex.h"
+
+#include "topology_mpi.h"
+#include "topology.h"
+
+std::vector<GModel*> PAR::createNewModels(GModel *gModel, GModel *global, int nparts, std::vector<int> myPart)
+{
+    int nbproc, myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &nbproc);
+    
+    int maxDim = -1;
+    std::vector<GModel*> newModels(myPart.size(), nullptr);
+    
+    for (unsigned int i = 0; i < myPart.size(); i++)
+    {
+        newModels[i] = new GModel();
+    }
+    
+    //Loop over regions
+    for(GModel::riter it = gModel->firstRegion(); it != gModel->lastRegion(); ++it)
+    {
+        GRegion *r = *it;
+        
+        std::vector<GRegion *> newModelHaveRegion(nparts, nullptr);
+        std::vector<GRegion *> globalHaveRegion(nparts, nullptr);
+        
+        //Tetrahedra
+        for(unsigned int i = 0; i < r->tetrahedra.size(); i++)
+        {
+            const int partition = r->tetrahedra[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+                
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->tetrahedra.push_back(r->tetrahedra[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->tetrahedra.push_back(r->tetrahedra[i]);
+        }
+        
+        
+        //Hexahedra
+        for(unsigned int i = 0; i < r->hexahedra.size(); i++)
+        {
+            const int partition = r->hexahedra[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+                
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->hexahedra.push_back(r->hexahedra[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->hexahedra.push_back(r->hexahedra[i]);
+        }
+        
+        //Prisms
+        for(unsigned int i = 0; i < r->prisms.size(); i++)
+        {
+            const int partition = r->prisms[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->prisms.push_back(r->prisms[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->prisms.push_back(r->prisms[i]);
+        }
+        
+        //Pyramids
+        for(unsigned int i = 0; i < r->pyramids.size(); i++)
+        {
+            const int partition = r->pyramids[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+                
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->pyramids.push_back(r->pyramids[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->pyramids.push_back(r->pyramids[i]);
+        }
+        
+        //Trihedra
+        for(unsigned int i = 0; i < r->trihedra.size(); i++)
+        {
+            const int partition = r->trihedra[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+                
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->trihedra.push_back(r->trihedra[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->trihedra.push_back(r->trihedra[i]);
+        }
+        
+        //Polyhedra
+        for(unsigned int i = 0; i < r->polyhedra.size(); i++)
+        {
+            const int partition = r->polyhedra[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveRegion[partition])
+                {
+                    discreteRegion *dr = new discreteRegion(newModels[position], newModels[position]->getNumRegions()+1);
+                    newModels[position]->add(dr);
+                    newModelHaveRegion[partition] = dr;
+                
+                    maxDim = 3;
+                    addPhysical(newModels[position], dr, gModel, r, partition, maxDim);
+                }
+            
+                newModelHaveRegion[partition]->polyhedra.push_back(r->polyhedra[i]);
+            }
+            if(!globalHaveRegion[partition] && global != nullptr)
+            {
+                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
+                global->add(drGlobal);
+                globalHaveRegion[partition] = drGlobal;
+                
+                maxDim = 3;
+                addPhysical(global, drGlobal, gModel, r, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveRegion[partition]->polyhedra.push_back(r->polyhedra[i]);
+        }
+    }
+    
+    //Loop over faces
+    for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
+    {
+        GFace *f = *it;
+        
+        std::vector<GFace *> newModelHaveFace(nparts, nullptr);
+        std::vector<GFace *> globalHaveFace(nparts, nullptr);
+        
+        //Triangles
+        for(unsigned int i = 0; i < f->triangles.size(); i++)
+        {
+            const int partition = f->triangles[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveFace[partition])
+                {
+                    discreteFace *df = new discreteFace(newModels[position], newModels[position]->getNumFaces()+1);
+                    newModels[position]->add(df);
+                    newModelHaveFace[partition] = df;
+                
+                    if(maxDim == -1)
+                    {
+                        maxDim = 2;
+                    }
+                    addPhysical(newModels[position], df, gModel, f, partition, maxDim);
+                }
+            
+                newModelHaveFace[partition]->triangles.push_back(f->triangles[i]);
+            }
+            if(!globalHaveFace[partition] && global != nullptr)
+            {
+                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
+                global->add(dfGlobal);
+                globalHaveFace[partition] = dfGlobal;
+                
+                if(maxDim == -1)
+                {
+                    maxDim = 2;
+                }
+                addPhysical(global, dfGlobal, gModel, f, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveFace[partition]->triangles.push_back(f->triangles[i]);
+        }
+        
+        //Quadrangles
+        for(unsigned int i = 0; i < f->quadrangles.size(); i++)
+        {
+            const int partition = f->quadrangles[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveFace[partition])
+                {
+                    discreteFace *df = new discreteFace(newModels[position], newModels[position]->getNumFaces()+1);
+                    newModels[position]->add(df);
+                    newModelHaveFace[partition] = df;
+                
+                    if(maxDim == -1)
+                    {
+                        maxDim = 2;
+                    }
+                    addPhysical(newModels[position], df, gModel, f, partition, maxDim);
+                }
+            
+                newModelHaveFace[partition]->quadrangles.push_back(f->quadrangles[i]);
+            }
+            if(!globalHaveFace[partition] && global != nullptr)
+            {
+                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
+                global->add(dfGlobal);
+                globalHaveFace[partition] = dfGlobal;
+                
+                if(maxDim == -1)
+                {
+                    maxDim = 2;
+                }
+                addPhysical(global, dfGlobal, gModel, f, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveFace[partition]->quadrangles.push_back(f->quadrangles[i]);
+        }
+        
+        //Polygons
+        for(unsigned int i = 0; i < f->polygons.size(); i++)
+        {
+            const int partition = f->polygons[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveFace[partition])
+                {
+                    discreteFace *df = new discreteFace(newModels[position], newModels[position]->getNumFaces()+1);
+                    newModels[position]->add(df);
+                    newModelHaveFace[partition] = df;
+                
+                    if(maxDim == -1)
+                    {
+                        maxDim = 2;
+                    }
+                    addPhysical(newModels[position], df, gModel, f, partition, maxDim);
+                }
+            
+                newModelHaveFace[partition]->polygons.push_back(f->polygons[i]);
+            }
+            if(!globalHaveFace[partition] && global != nullptr)
+            {
+                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
+                global->add(dfGlobal);
+                globalHaveFace[partition] = dfGlobal;
+                
+                if(maxDim == -1)
+                {
+                    maxDim = 2;
+                }
+                addPhysical(global, dfGlobal, gModel, f, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveFace[partition]->polygons.push_back(f->polygons[i]);
+        }
+    }
+    
+    //Loop over edges
+    for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
+    {
+        GEdge *e = *it;
+        
+        std::vector<GEdge *> newModelHaveEdge(nparts, nullptr);
+        std::vector<GEdge *> globalHaveEdge(nparts, nullptr);
+        
+        //Lines
+        for(unsigned int i = 0; i < e->lines.size(); i++)
+        {
+            const int partition = e->lines[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveEdge[partition])
+                {
+                    discreteEdge *de = new discreteEdge(newModels[position], newModels[position]->getNumEdges()+1, nullptr, nullptr);
+                    newModels[position]->add(de);
+                    newModelHaveEdge[partition] = de;
+                    
+                    if(maxDim == -1)
+                    {
+                        maxDim = 1;
+                    }
+                    addPhysical(newModels[position], de, gModel, e, partition, maxDim);
+                }
+            
+                newModelHaveEdge[partition]->lines.push_back(e->lines[i]);
+            }
+            if(!globalHaveEdge[partition] && global != nullptr)
+            {
+                discreteEdge *deGlobal = new discreteEdge(global, global->getNumEdges()+1, nullptr, nullptr);
+                global->add(deGlobal);
+                globalHaveEdge[partition] = deGlobal;
+                
+                if(maxDim == -1)
+                {
+                    maxDim = 1;
+                }
+                addPhysical(global, deGlobal, gModel, e, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveEdge[partition]->lines.push_back(e->lines[i]);
+        }
+    }
+    
+    //Loop over vertices
+    for(GModel::viter it = gModel->firstVertex(); it != gModel->lastVertex(); ++it)
+    {
+        GVertex *v = *it;
+        
+        std::vector<GVertex *> newModelHaveVertex(nparts, nullptr);
+        std::vector<GVertex *> globalHaveVertex(nparts, nullptr);
+        
+        for(unsigned int i = 0; i < v->points.size(); i++)
+        {
+            const int partition = v->points[i]->getPartition();
+            if(std::find(myPart.begin(), myPart.end(), partition) != myPart.end())
+            {
+                int position = getPos(myPart, partition);
+                if(!newModelHaveVertex[partition])
+                {
+                    discreteVertex *dv = new discreteVertex(newModels[position], newModels[position]->getNumVertices()+1);
+                    newModels[position]->add(dv);
+                    newModelHaveVertex[partition] = dv;
+                
+                    if(maxDim == -1)
+                    {
+                        maxDim = 0;
+                    }
+                    addPhysical(newModels[position], dv, gModel, v, partition, maxDim);
+                }
+            
+                newModelHaveVertex[partition]->points.push_back(v->points[i]);
+            }
+            if(!globalHaveVertex[partition] && global != nullptr)
+            {
+                discreteVertex *dvGlobal = new discreteVertex(global, global->getNumVertices()+1);
+                global->add(dvGlobal);
+                globalHaveVertex[partition] = dvGlobal;
+                
+                if(maxDim == -1)
+                {
+                    maxDim = 0;
+                }
+                addPhysical(global, dvGlobal, gModel, v, partition, maxDim);
+            }
+            
+            if(global != nullptr) globalHaveVertex[partition]->points.push_back(v->points[i]);
+        }
+    }
+    
+    return newModels;
+}
+
+void PAR::addPhysical(GModel *newModel, GEntity *newEntity, GModel *oldModel, GEntity *oldEntity, int partition, int maxDim)
+{
+    std::vector<int> oldPhysical = oldEntity->getPhysicalEntities();
+    std::string name;
+    
+    if(maxDim == oldEntity->dim())
+    {
+        name = "_omega{";
+        name += std::to_string(partition);
+        name += "}";
+        
+        const int number = oldModel->getMaxPhysicalNumber(-1)+1+partition;
+        
+        newModel->setPhysicalName(name, newEntity->dim(), number);
+        newEntity->addPhysicalEntity(number);
+    }
+    
+    for(unsigned int i = 0; i < oldPhysical.size(); i++)
+    {
+        name = oldModel->getPhysicalName(oldEntity->dim(), oldPhysical[i]);
+        
+        if(name[0] != '_')
+        {
+            newModel->setPhysicalName(name, newEntity->dim(), oldPhysical[i]);
+            newEntity->addPhysicalEntity(oldPhysical[i]);
+        }
+    }
+}
+
+int PAR::getPos(std::vector<int> myPart, int partition)
+{
+    for(unsigned int i = 0; i < myPart.size(); i++)
+    {
+        if(myPart[i] == partition) return i;
+    }
+    
+    return -1;
+}

Partitioning/topology_mpi.h

+#ifndef TOPOLOGY_MPI_INCLUDED
+#define TOPOLOGY_MPI_INCLUDED
+
+#include <unordered_map>
+
+#include "GModel.h"
+
+#include "partitionFace.h"
+#include "partitionEdge.h"
+#include "partitionVertex.h"
+
+#include "MFaceHash.h"
+#include "MEdgeHash.h"
+
+namespace PAR {
+    std::vector<GModel*> createNewModels(GModel *gModel, GModel *global, int nparts, std::vector<int> myPart);
+    void addPhysical(GModel *newModel, GEntity *newEntity, GModel *oldModel, GEntity *oldEntity, int partition, int maxDim);
+    int getPos(std::vector<int> myPart, int partition);
+}
+
+#endif //TOPOLOGY_MPI_INCLUDED