diff --git a/Partitioning/CMakeLists.txt b/Partitioning/CMakeLists.txt
deleted file mode 100755
index d69d88ea4334f2040a66cc6dbdd3d53567969e82..0000000000000000000000000000000000000000
--- a/Partitioning/CMakeLists.txt
+++ /dev/null
@@ -1,51 +0,0 @@
-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)
diff --git a/Partitioning/graph.cpp b/Partitioning/graph.cpp
deleted file mode 100755
index 30432292e6d7db5a7688cb0f266b0f01ddfe4a4c..0000000000000000000000000000000000000000
--- a/Partitioning/graph.cpp
+++ /dev/null
@@ -1,160 +0,0 @@
-#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
-
diff --git a/Partitioning/graph.h b/Partitioning/graph.h
deleted file mode 100755
index d57706f7c225c3b39cbc4a011f2adb2c5aabff2b..0000000000000000000000000000000000000000
--- a/Partitioning/graph.h
+++ /dev/null
@@ -1,19 +0,0 @@
-#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
diff --git a/Partitioning/io.cpp b/Partitioning/io.cpp
deleted file mode 100644
index 7021a881c2ddd525fd332b865512944ae77a17a7..0000000000000000000000000000000000000000
--- a/Partitioning/io.cpp
+++ /dev/null
@@ -1,374 +0,0 @@
-#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;
-}
-
-
-
diff --git a/Partitioning/io.h b/Partitioning/io.h
deleted file mode 100644
index 1dfa6c9319e055cf650a3e550abc917f2183070d..0000000000000000000000000000000000000000
--- a/Partitioning/io.h
+++ /dev/null
@@ -1,11 +0,0 @@
-#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 */
diff --git a/Partitioning/main.cpp b/Partitioning/main.cpp
deleted file mode 100755
index 3c943c997c7f3f0438652866f3e0a9a86a03a6b2..0000000000000000000000000000000000000000
--- a/Partitioning/main.cpp
+++ /dev/null
@@ -1,233 +0,0 @@
-#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;
-}
diff --git a/Partitioning/main_mpi.cpp b/Partitioning/main_mpi.cpp
deleted file mode 100644
index c3b75dd73e30e48b18e4504c530f732623a0105e..0000000000000000000000000000000000000000
--- a/Partitioning/main_mpi.cpp
+++ /dev/null
@@ -1,200 +0,0 @@
-#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;
-}
diff --git a/Partitioning/main_mpi.h b/Partitioning/main_mpi.h
deleted file mode 100644
index 15da23ac2dffb59c740e7022da527d211b8ed88f..0000000000000000000000000000000000000000
--- a/Partitioning/main_mpi.h
+++ /dev/null
@@ -1,8 +0,0 @@
-#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
diff --git a/Partitioning/topology.cpp b/Partitioning/topology.cpp
deleted file mode 100644
index 656d42f93921ba5ca5e28c1cdfa99c5fcf0fd48c..0000000000000000000000000000000000000000
--- a/Partitioning/topology.cpp
+++ /dev/null
@@ -1,1190 +0,0 @@
-#include <iostream>
-#include <set>
-#include <map>
-#include <unordered_map>
-#include <string>
-
-#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.h"
-
-int SEQ::createPartitionBoundaries(GModel *model, std::vector<GModel*> models)
-{
-    unsigned int numElem[6];
-    const int meshDim = model->getNumMeshElements(numElem);
-    
-    std::set<partitionFace*, Less_partitionFace> pfaces;
-    std::set<partitionEdge*, Less_partitionEdge> pedges;
-    std::set<partitionVertex*, Less_partitionVertex> pvertices;
-    
-    std::vector< std::set<partitionFace*, Less_partitionFace> > pfacesOfModels(models.size(), std::set<partitionFace*, Less_partitionFace>() );
-    std::vector< std::set<partitionEdge*, Less_partitionEdge> > pedgesOfModels(models.size(), std::set<partitionEdge*, Less_partitionEdge>() );
-    std::vector< std::set<partitionVertex*, Less_partitionVertex> > pverticesOfModels(models.size(), std::set<partitionVertex*, Less_partitionVertex>() );
-
-    std::unordered_map<MFace, std::vector<MElement*> , Hash_Face, Equal_Face> faceToElement;
-    std::unordered_map<MEdge, std::vector<MElement*> , Hash_Edge, Equal_Edge> edgeToElement;
-    std::unordered_map<MVertex*, std::vector<MElement*> > vertexToElement;
-    
-    //Create partition faces
-    std::cout << "\tCreate partition faces... " << std::flush;
-    if (meshDim == 3)
-    {
-        for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it)
-        {
-            fillit_(faceToElement, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
-            fillit_(faceToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
-            fillit_(faceToElement, (*it)->prisms.begin(), (*it)->prisms.end());
-            fillit_(faceToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
-            fillit_(faceToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
-            fillit_(faceToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
-        }
-        
-        for(std::unordered_map<MFace, std::vector<MElement*> , Hash_Face, Equal_Face>::const_iterator it = faceToElement.begin(); it != faceToElement.end(); ++it)
-        {
-            MFace f = it->first;
-            std::vector<MElement*> voe = it->second;
-            
-            assignPartitionBoundary(model, f, pfaces, voe, models, pfacesOfModels);
-        }
-    }
-    std::cout << "Done!" << std::endl;
-    
-    //Create partition edges
-    std::cout << "\tCreate partition edges... " << std::flush;
-    if (meshDim >= 2)
-    {
-        if (meshDim == 2)
-        {
-            for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it)
-            {
-                fillit_(edgeToElement, (*it)->triangles.begin(), (*it)->triangles.end());
-                fillit_(edgeToElement, (*it)->quadrangles.begin(), (*it)->quadrangles.end());
-                fillit_(edgeToElement, (*it)->polygons.begin(), (*it)->polygons.end());
-            }
-        }
-        
-        if (meshDim == 3)
-        {
-            for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it)
-            {
-                fillit_(edgeToElement, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
-                fillit_(edgeToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
-                fillit_(edgeToElement, (*it)->prisms.begin(), (*it)->prisms.end());
-                fillit_(edgeToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
-                fillit_(edgeToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
-                fillit_(edgeToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
-            }
-        }
-        
-        for(std::unordered_map<MEdge, std::vector<MElement*> , Hash_Edge, Equal_Edge>::const_iterator it = edgeToElement.begin(); it != edgeToElement.end(); ++it)
-        {
-            MEdge e = it->first;
-            
-            std::vector<MElement*> voe = it->second;
-            
-            assignPartitionBoundary(model, e, pedges, voe, pfaces, models, pedgesOfModels);
-        }
-    }
-    std::cout << "Done!" << std::endl;
-    
-    //Create partition vertices
-    std::cout << "\tCreate partition vertices... "  << std::flush;
-    if (meshDim >= 1)
-    {
-        if (meshDim == 1)
-        {
-            for(GModel::eiter it = model->firstEdge(); it != model->lastEdge(); ++it)
-            {
-                fillit_(vertexToElement, (*it)->lines.begin(), (*it)->lines.end());
-            }
-        }
-        
-        if (meshDim == 2)
-        {
-            for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it)
-            {
-                fillit_(vertexToElement, (*it)->triangles.begin(), (*it)->triangles.end());
-                fillit_(vertexToElement, (*it)->quadrangles.begin(), (*it)->quadrangles.end());
-                fillit_(vertexToElement, (*it)->polygons.begin(), (*it)->polygons.end());
-            }
-        }
-        
-        if (meshDim == 3)
-        {
-            for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it)
-            {
-                fillit_(vertexToElement, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
-                fillit_(vertexToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
-                fillit_(vertexToElement, (*it)->prisms.begin(), (*it)->prisms.end());
-                fillit_(vertexToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
-                fillit_(vertexToElement, (*it)->trihedra.begin(), (*it)->trihedra.end());
-                fillit_(vertexToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
-            }
-        }
-        
-        for(std::unordered_map<MVertex*, std::vector<MElement*> >::const_iterator it = vertexToElement.begin(); it != vertexToElement.end(); ++it)
-        {
-            MVertex *v = it->first;
-            std::vector<MElement*> voe = it->second;
-            
-            assignPartitionBoundary(model, v, pvertices, voe, pedges, pfaces, models, pverticesOfModels);
-        }
-    }
-    std::cout << "Done!" << std::endl;
-    
-    return 0;
-}
-
-template <class ITERATOR>
-void SEQ::fillit_(std::unordered_map<MFace, std::vector<MElement*> , Hash_Face, Equal_Face> &faceToElement, ITERATOR it_beg, ITERATOR it_end)
-{
-    for (ITERATOR IT = it_beg; IT != it_end ; ++IT)
-    {
-        MElement *el = *IT;
-        for(unsigned int j = 0; j < el->getNumFaces(); j++)
-        {
-            faceToElement[el->getFace(j)].push_back(el);
-        }
-    }
-}
-
-template <class ITERATOR>
-void SEQ::fillit_(std::unordered_map<MEdge, std::vector<MElement*> , Hash_Edge, Equal_Edge> &edgeToElement, ITERATOR it_beg, ITERATOR it_end)
-{
-    for (ITERATOR IT = it_beg; IT != it_end; ++IT)
-    {
-        MElement *el = *IT;
-        for(unsigned int j = 0; j < el->getNumEdges(); j++)
-        {
-            edgeToElement[el->getEdge(j)].push_back(el);
-        }
-    }
-}
-
-template <class ITERATOR>
-void SEQ::fillit_(std::unordered_map<MVertex*, std::vector<MElement*> > &vertexToElement, ITERATOR it_beg, ITERATOR it_end)
-{
-    for (ITERATOR IT = it_beg; IT != it_end ; ++IT)
-    {
-        MElement *el = *IT;
-        for(unsigned int j = 0; j < el->getNumVertices(); j++)
-        {
-            vertexToElement[el->getVertex(j)].push_back(el);
-        }
-    }
-}
-
-void SEQ::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)
-{
-    std::vector<int> v2;
-    v2.push_back(v[0]->getPartition());
-    
-    for (unsigned int i = 1; i < v.size(); i++)
-    {
-        bool found = false;
-        for (unsigned int j = 0; j < v2.size(); j++)
-        {
-            if (v[i]->getPartition() == v2[j])
-            {
-                found = true;
-                break;
-            }
-        }
-        
-        if (!found)
-        {
-            v2.push_back(v[i]->getPartition());
-        }
-    }
-    
-    //If there is less than two partitions touching the face: stop
-    if (v2.size() < 2)
-    {
-        return;
-    }
-    
-    const int numPhysical = model->getMaxPhysicalNumber(-1)+1;
-    
-    partitionFace pf(model, 1, v2);
-    std::set<partitionFace*, Less_partitionFace>::iterator it = pfaces.find(&pf);
-    
-    partitionFace *ppf;
-     //Create the new partition entity for the global mesh
-    if (it == pfaces.end())
-    {
-        //Create new entity and add them to the global model
-        ppf = new  partitionFace(model, -(int)pfaces.size()-1, v2);
-        pfaces.insert(ppf);
-        model->add(ppf);
-        
-        //Create its new physical name
-        ppf->addPhysicalEntity(numPhysical);
-        
-        std::string name = "_sigma{";
-        for(unsigned int j = 0; j < ppf->_partitions.size(); j++)
-        {
-            if(j > 0)
-            {
-                name += ",";
-            }
-            name += std::to_string(ppf->_partitions[j]);
-        }
-        name += "}";
-        
-        model->setPhysicalName(name, ppf->dim(), numPhysical);
-    }
-    else
-    {
-        ppf = *it;
-    }
-    
-    int numFace = 0;
-    for(unsigned int i = 0; i < v[0]->getNumFaces(); i++)
-    {
-        const MFace e = v[0]->getFace(i);
-        if(e == me)
-        {
-            numFace = i;
-            break;
-        }
-    }
-    
-    if(me.getNumVertices() == 3)
-    {
-        std::vector<MVertex*> verts;
-        v[0]->getFaceVertices(numFace, verts);
-        
-        if(verts.size() == 3)
-        {
-            ppf->triangles.push_back(new MTriangle(verts));
-        }
-        else if(verts.size() == 6)
-        {
-            ppf->triangles.push_back(new MTriangle6(verts));
-        }
-        else
-        {
-            ppf->triangles.push_back(new MTriangleN(verts, verts[0]->getPolynomialOrder()));
-        }
-    }
-    else if(me.getNumVertices() == 4)
-    {
-        std::vector<MVertex*> verts;
-        v[0]->getFaceVertices(numFace, verts);
-        
-        if(verts.size() == 4)
-        {
-            ppf->quadrangles.push_back(new MQuadrangle(verts));
-        }
-        else if(verts.size() == 8)
-        {
-            ppf->quadrangles.push_back(new MQuadrangle8(verts));
-        }
-        else if(verts.size() == 9)
-        {
-            ppf->quadrangles.push_back(new MQuadrangle9(verts));
-        }
-        else
-        {
-            ppf->quadrangles.push_back(new MQuadrangleN(verts, verts[0]->getPolynomialOrder()));
-        }
-    }
-    
-    for(unsigned int i = 0; i < v2.size(); i++)
-    {
-        partitionFace pf(model, 1, v2);
-        std::set<partitionFace*, Less_partitionFace>::iterator it = pfacesOfModels[v2[i]].find(&pf);
-        
-        //Create the new partition entity for partitioned meshes
-        partitionFace *ppf;
-        if (it == pfacesOfModels[v2[i]].end())
-        {
-            //Create new entity and add them to partitioned models
-            ppf = new  partitionFace(models[v2[i]], -(int)pfacesOfModels[v2[i]].size()-1, v2);
-            pfacesOfModels[v2[i]].insert(ppf);
-            models[v2[i]]->add(ppf);
-            
-            //Create its new physical name
-            ppf->addPhysicalEntity(numPhysical);
-            
-            std::string name = "_sigma{";
-            for(unsigned int j = 0; j < ppf->_partitions.size(); j++)
-            {
-                if(j > 0)
-                {
-                    name += ",";
-                }
-                name += std::to_string(ppf->_partitions[j]);
-            }
-            name += "}";
-            
-            models[v2[i]]->setPhysicalName(name, ppf->dim(), numPhysical);
-        }
-        else
-        {
-            ppf = *it;
-        }
-        
-        int numEdge = 0;
-        int numElm = 0;
-        for(unsigned int j = 0; j < v.size(); j++)
-        {
-            if(v[j]->getPartition() == v2[i])
-            {
-                numElm = j;
-                for(unsigned int k = 0; k < v[j]->getNumFaces(); k++)
-                {
-                    const MFace e = v[j]->getFace(k);
-                    if(e == me)
-                    {
-                        numFace = k;
-                        break;
-                    }
-                }
-                break;
-            }
-        }
-        
-        if(me.getNumVertices() == 3)
-        {
-            std::vector<MVertex*> verts;
-            v[numElm]->getFaceVertices(numFace, verts);
-            
-            if(verts.size() == 3)
-            {
-                ppf->triangles.push_back(new MTriangle(verts));
-            }
-            else if(verts.size() == 6)
-            {
-                ppf->triangles.push_back(new MTriangle6(verts));
-            }
-            else
-            {
-                ppf->triangles.push_back(new MTriangleN(verts, verts[0]->getPolynomialOrder()));
-            }
-        }
-        else if(me.getNumVertices() == 4)
-        {
-            std::vector<MVertex*> verts;
-            v[numElm]->getFaceVertices(numFace, verts);
-            
-            if(verts.size() == 4)
-            {
-                ppf->quadrangles.push_back(new MQuadrangle(verts));
-            }
-            else if(verts.size() == 8)
-            {
-                ppf->quadrangles.push_back(new MQuadrangle8(verts));
-            }
-            else if(verts.size() == 9)
-            {
-                ppf->quadrangles.push_back(new MQuadrangle9(verts));
-            }
-            else
-            {
-                ppf->quadrangles.push_back(new MQuadrangleN(verts, verts[0]->getPolynomialOrder()));
-            }
-        }
-    }
-}
-
-void SEQ::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)
-{
-    std::vector<int> v2;
-    v2.push_back(v[0]->getPartition());
-    
-    for (unsigned int i = 1; i < v.size(); i++)
-    {
-        bool found = false;
-        for (unsigned int j = 0; j < v2.size(); j++)
-        {
-            if (v[i]->getPartition() == v2[j])
-            {
-                found = true;
-                break;
-            }
-        }
-        
-        if (!found)
-        {
-            v2.push_back(v[i]->getPartition());
-        }
-    }
-    
-    //If there is less than two partitions touching the edge: stop
-    if (v2.size() < 2)
-    {
-        return;
-    }
-    
-    partitionFace pf(model, 1, v2);
-    std::set<partitionFace*, Less_partitionFace>::iterator itf = pfaces.find(&pf);
-    
-    //If the edge is on a partitionFace
-    if (itf != pfaces.end())
-    {
-        return;
-    }
-    
-    const int numPhysical = model->getMaxPhysicalNumber(-1)+1;
-    
-    partitionEdge pe(model, 1, nullptr, nullptr, v2);
-    std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedges.find(&pe);
-    
-    partitionEdge *ppe;
-    //Create the new partition entity for the global mesh
-    if (it == pedges.end())
-    {
-        //Create new entity and add them to the global model
-        ppe = new  partitionEdge(model, -(int)pedges.size()-1, 0, 0, v2);
-        pedges.insert(ppe);
-        model->add(ppe);
-        
-        //Create its new physical name
-        ppe->addPhysicalEntity(numPhysical);
-        
-        std::string name = "_sigma{";
-        for(unsigned int j = 0; j < ppe->_partitions.size(); j++)
-        {
-            if(j > 0)
-            {
-                name += ",";
-            }
-            name += std::to_string(ppe->_partitions[j]);
-        }
-        name += "}";
-        
-        model->setPhysicalName(name, ppe->dim(), numPhysical);
-    }
-    else
-    {
-        ppe = *it;
-    }
-    
-    int numEdge = 0;
-    for(unsigned int i = 0; i < v[0]->getNumEdges(); i++)
-    {
-        const MEdge e = v[0]->getEdge(i);
-        if(e == me)
-        {
-            numEdge = i;
-            break;
-        }
-    }
-    
-    if(me.getNumVertices() == 2)
-    {
-        std::vector<MVertex*> verts;
-        v[0]->getEdgeVertices(numEdge, verts);
-        
-        if(verts.size() == 2)
-        {
-            ppe->lines.push_back(new MLine(verts));
-        }
-        else if(verts.size() == 3)
-        {
-            ppe->lines.push_back(new MLine3(verts));
-        }
-        else
-        {
-            ppe->lines.push_back(new MLineN(verts));
-        }
-    }
-    
-    for(unsigned int i = 0; i < v2.size(); i++)
-    {
-        partitionEdge pe(models[v2[i]], 1, nullptr, nullptr, v2);
-        std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedgesOfModels[v2[i]].find(&pe);
-        
-        //Create the new partition entity for partitioned meshes
-        partitionEdge *ppe;
-        if (it == pedgesOfModels[v2[i]].end())
-        {
-            //Create new entity and add them to partitioned models
-            ppe = new  partitionEdge(models[v2[i]], -(int)pedgesOfModels[v2[i]].size()-1, 0, 0, v2);
-            pedgesOfModels[v2[i]].insert(ppe);
-            models[v2[i]]->add(ppe);
-            
-            //Create its new physical name
-            ppe->addPhysicalEntity(numPhysical);
-            
-            std::string name = "_sigma{";
-            for(unsigned int j = 0; j < ppe->_partitions.size(); j++)
-            {
-                if(j > 0)
-                {
-                    name += ",";
-                }
-                name += std::to_string(ppe->_partitions[j]);
-            }
-            name += "}";
-            
-            models[v2[i]]->setPhysicalName(name, ppe->dim(), numPhysical);
-        }
-        else
-        {
-            ppe = *it;
-        }
-        
-        int numEdge = 0;
-        int numElm = 0;
-
-        for(unsigned int j = 0; j < v.size(); j++)
-        {
-            if(v[j]->getPartition() == v2[i])
-            {
-                numElm = j;
-                for(unsigned int k = 0; k < v[j]->getNumEdges(); k++)
-                {
-                    const MEdge e = v[j]->getEdge(k);
-                    if(e == me)
-                    {
-                        numEdge = k;
-                        break;
-                    }
-                }
-                break;
-            }
-        }
-        
-        if(me.getNumVertices() == 2)
-        {
-            std::vector<MVertex*> verts;
-            v[numElm]->getEdgeVertices(numEdge, verts);
-            
-            if(verts.size() == 2)
-            {
-                ppe->lines.push_back(new MLine(verts));
-            }
-            else if(verts.size() == 3)
-            {
-                ppe->lines.push_back(new MLine3(verts));
-            }
-            else
-            {
-                ppe->lines.push_back(new MLineN(verts));
-            }
-        }
-    }
-}
-
-void SEQ::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<int> v2;
-    v2.push_back(v[0]->getPartition());
-    
-    for (unsigned int i = 1; i < v.size(); i++)
-    {
-        bool found = false;
-        for (unsigned int j = 0; j < v2.size(); j++)
-        {
-            if (v[i]->getPartition() == v2[j])
-            {
-                found = true;
-                break;
-            }
-        }
-        
-        if (!found)
-        {
-            v2.push_back(v[i]->getPartition());
-        }
-    }
-    
-    //If there is less than two partitions touching the vertex: stop
-    if (v2.size() < 2)
-    {
-        return;
-    }
-    
-    partitionFace pf(model, 1, v2);
-    std::set<partitionFace*, Less_partitionFace>::iterator itf = pfaces.find(&pf);
-    
-    //If the vertex is on a partitionFace
-    if (itf != pfaces.end())
-    {
-        return;
-    }
-    
-    partitionEdge pe(model, 1, 0, 0, v2);
-    std::set<partitionEdge*, Less_partitionEdge>::iterator ite = pedges.find(&pe);
-    
-    //If the vertex is on a partitionFace
-    if (ite != pedges.end())
-    {
-        return;
-    }
-
-    const int numPhysical = model->getMaxPhysicalNumber(-1)+1;
-    
-    partitionVertex pv(model, 1, v2);
-    std::set<partitionVertex*, Less_partitionVertex>::iterator it = pvertices.find(&pv);
-    
-    partitionVertex *ppv;
-    //Create the new partition entity for the global mesh
-    if (it == pvertices.end())
-    {
-        ppv = new partitionVertex(model, -(int)pvertices.size()-1,v2);
-        pvertices.insert (ppv);
-        model->add(ppv);
-        
-        //Create its new physical name
-        ppv->addPhysicalEntity(numPhysical);
-        
-        std::string name = "_sigma{";
-        for(unsigned int j = 0; j < ppv->_partitions.size(); j++)
-        {
-            if(j > 0)
-            {
-                name += ",";
-            }
-            name += std::to_string(ppv->_partitions[j]);
-        }
-        name += "}";
-        
-        model->setPhysicalName(name, ppv->dim(), numPhysical);
-    }
-    else
-    {
-        ppv = *it;
-    }
-    
-    ppv->points.push_back(new MPoint(ve));
-    
-    for(unsigned int i = 0; i < v2.size(); i++)
-    {
-        partitionVertex pv(model, 1, v2);
-        std::set<partitionVertex*, Less_partitionVertex>::iterator it = pverticesOfModels[v2[i]].find(&pv);
-        
-        partitionVertex *ppv;
-        //Create the new partition entity for the global mesh
-        if (it == pverticesOfModels[v2[i]].end())
-        {
-            ppv = new partitionVertex(model, -(int)pverticesOfModels[v2[i]].size()-1,v2);
-            pverticesOfModels[v2[i]].insert (ppv);
-            models[v2[i]]->add(ppv);
-            
-            //Create its new physical name
-            ppv->addPhysicalEntity(numPhysical);
-            
-            std::string name = "_sigma{";
-            for(unsigned int j = 0; j < ppv->_partitions.size(); j++)
-            {
-                if(j > 0)
-                {
-                    name += ",";
-                }
-                name += std::to_string(ppv->_partitions[j]);
-            }
-            name += "}";
-            
-            models[v2[i]]->setPhysicalName(name, ppv->dim(), numPhysical);
-        }
-        else
-        {
-            ppv = *it;
-        }
-        
-        ppv->points.push_back(new MPoint(ve));
-    }
-}
-
-std::vector<GModel*> SEQ::createNewModels(GModel *gModel, GModel *global, int nparts)
-{
-    int maxDim = -1;
-    std::vector<GModel*> newModels(nparts, nullptr);
-    
-    for (unsigned int i = 0; i < nparts; 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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->tetrahedra.push_back(r->tetrahedra[i]);
-            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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->hexahedra.push_back(r->hexahedra[i]);
-            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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->prisms.push_back(r->prisms[i]);
-            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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->pyramids.push_back(r->pyramids[i]);
-            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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->trihedra.push_back(r->trihedra[i]);
-            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(!newModelHaveRegion[partition])
-            {
-                discreteRegion *dr = new discreteRegion(newModels[partition], newModels[partition]->getNumRegions()+1);
-                newModels[partition]->add(dr);
-                newModelHaveRegion[partition] = dr;
-                
-                discreteRegion *drGlobal = new discreteRegion(global, global->getNumRegions()+1);
-                global->add(drGlobal);
-                globalHaveRegion[partition] = drGlobal;
-                
-                maxDim = 3;
-                addPhysical(newModels[partition], dr, gModel, r, global, drGlobal, partition, maxDim);
-            }
-            
-            newModelHaveRegion[partition]->polyhedra.push_back(r->polyhedra[i]);
-            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 *> globalModelHaveFace(nparts, nullptr);
-        
-        //Triangles
-        for(unsigned int i = 0; i < f->triangles.size(); i++)
-        {
-            const int partition = f->triangles[i]->getPartition();
-            if(!newModelHaveFace[partition])
-            {
-                discreteFace *df = new discreteFace(newModels[partition], newModels[partition]->getNumFaces()+1);
-                newModels[partition]->add(df);
-                newModelHaveFace[partition] = df;
-                
-                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
-                global->add(dfGlobal);
-                globalModelHaveFace[partition] = dfGlobal;
-                
-                if(maxDim == -1)
-                {
-                    maxDim = 2;
-                }
-                addPhysical(newModels[partition], df, gModel, f, global, dfGlobal, partition, maxDim);
-            }
-            
-            newModelHaveFace[partition]->triangles.push_back(f->triangles[i]);
-            globalModelHaveFace[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(!newModelHaveFace[partition])
-            {
-                discreteFace *df = new discreteFace(newModels[partition], newModels[partition]->getNumFaces()+1);
-                newModels[partition]->add(df);
-                newModelHaveFace[partition] = df;
-                
-                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
-                global->add(dfGlobal);
-                globalModelHaveFace[partition] = dfGlobal;
-                
-                if(maxDim == -1)
-                {
-                    maxDim = 2;
-                }
-                addPhysical(newModels[partition], df, gModel, f, global, dfGlobal, partition, maxDim);
-            }
-            
-            newModelHaveFace[partition]->quadrangles.push_back(f->quadrangles[i]);
-            globalModelHaveFace[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(!newModelHaveFace[partition])
-            {
-                discreteFace *df = new discreteFace(newModels[partition], newModels[partition]->getNumFaces()+1);
-                newModels[partition]->add(df);
-                newModelHaveFace[partition] = df;
-                
-                discreteFace *dfGlobal = new discreteFace(global, global->getNumFaces()+1);
-                global->add(dfGlobal);
-                globalModelHaveFace[partition] = dfGlobal;
-                
-                if(maxDim == -1)
-                {
-                    maxDim = 2;
-                }
-                addPhysical(newModels[partition], df, gModel, f, global, dfGlobal, partition, maxDim);
-            }
-            
-            newModelHaveFace[partition]->polygons.push_back(f->polygons[i]);
-            globalModelHaveFace[partition]->quadrangles.push_back(f->quadrangles[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 *> globalModelHaveEdge(nparts, nullptr);
-        
-        //Lines
-        for(unsigned int i = 0; i < e->lines.size(); i++)
-        {
-            const int partition = e->lines[i]->getPartition();
-            if(!newModelHaveEdge[partition])
-            {
-                discreteEdge *de = new discreteEdge(newModels[partition], newModels[partition]->getNumEdges()+1, nullptr, nullptr);
-                newModels[partition]->add(de);
-                newModelHaveEdge[partition] = de;
-                
-                discreteEdge *deGlobal = new discreteEdge(global, global->getNumEdges()+1, nullptr, nullptr);
-                global->add(deGlobal);
-                globalModelHaveEdge[partition] = deGlobal;
-                
-                if(maxDim == -1)
-                {
-                    maxDim = 1;
-                }
-                addPhysical(newModels[partition], de, gModel, e, global, deGlobal, partition, maxDim);
-            }
-            
-            newModelHaveEdge[partition]->lines.push_back(e->lines[i]);
-            globalModelHaveEdge[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 *> globalModelHaveVertex(nparts, nullptr);
-        
-        for(unsigned int i = 0; i < v->points.size(); i++)
-        {
-            const int partition = v->points[i]->getPartition();
-            if(!newModelHaveVertex[partition])
-            {
-                discreteVertex *dv = new discreteVertex(newModels[partition], newModels[partition]->getNumVertices()+1);
-                newModels[partition]->add(dv);
-                newModelHaveVertex[partition] = dv;
-                
-                discreteVertex *dvGlobal = new discreteVertex(global, global->getNumVertices()+1);
-                global->add(dvGlobal);
-                globalModelHaveVertex[partition] = dvGlobal;
-                
-                if(maxDim == -1)
-                {
-                    maxDim = 0;
-                }
-                addPhysical(newModels[partition], dv, gModel, v, global, dvGlobal, partition, maxDim);
-            }
-            
-            newModelHaveVertex[partition]->points.push_back(v->points[i]);
-            globalModelHaveVertex[partition]->points.push_back(v->points[i]);
-        }
-    }
-    
-    return newModels;
-}
-
-void SEQ::assignMeshVerticesToModel(GModel *gModel)
-{
-    std::unordered_map<MVertex*, GEntity*> vertexToEntity;
-    
-    //Loop over regions
-    for(GModel::riter it = gModel->firstRegion(); it != gModel->lastRegion(); ++it)
-    {
-        GRegion *r = *it;
-        
-        fillVertexToEntity(vertexToEntity, r, r->tetrahedra.begin(), r->tetrahedra.end());
-        fillVertexToEntity(vertexToEntity, r, r->hexahedra.begin(), r->hexahedra.end());
-        fillVertexToEntity(vertexToEntity, r, r->prisms.begin(), r->prisms.end());
-        fillVertexToEntity(vertexToEntity, r, r->pyramids.begin(), r->pyramids.end());
-        fillVertexToEntity(vertexToEntity, r, r->trihedra.begin(), r->trihedra.end());
-        fillVertexToEntity(vertexToEntity, r, r->polyhedra.begin(), r->polyhedra.end());
-    }
-    
-    //Loop over faces
-    for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
-    {
-        GFace *f = *it;
-        
-        fillVertexToEntity(vertexToEntity, f, f->triangles.begin(), f->triangles.end());
-        fillVertexToEntity(vertexToEntity, f, f->quadrangles.begin(), f->quadrangles.end());
-        fillVertexToEntity(vertexToEntity, f, f->polygons.begin(), f->polygons.end());
-    }
-    
-    //Loop over edges
-    for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
-    {
-        GEdge *e = *it;
-        
-        fillVertexToEntity(vertexToEntity, e, e->lines.begin(), e->lines.end());
-    }
-    
-    //Loop over vertices
-    for(GModel::viter it = gModel->firstVertex(); it != gModel->lastVertex(); ++it)
-    {
-        GVertex *v = *it;
-        
-        fillVertexToEntity(vertexToEntity, v, v->points.begin(), v->points.end());
-    }
-    
-    //Fill the entities
-    for(std::unordered_map<MVertex*, GEntity*>::iterator it = vertexToEntity.begin(); it != vertexToEntity.end(); ++it)
-    {
-        it->second->addMeshVertex(it->first);
-    }
-}
-
-template <class ITERATOR>
-void SEQ::fillVertexToEntity(std::unordered_map<MVertex*, GEntity*> &vertexToEntity, GEntity* entity, ITERATOR it_beg, ITERATOR it_end)
-{
-    for(ITERATOR it = it_beg; it != it_end; ++it)
-    {
-        const int numVertices = (*it)->getNumVertices();
-        for(unsigned int j = 0; j < numVertices; j++)
-        {
-            if(vertexToEntity.find((*it)->getVertex(j)) != vertexToEntity.end())
-            {
-                if(vertexToEntity[(*it)->getVertex(j)]->dim() > entity->dim())
-                {
-                    vertexToEntity[(*it)->getVertex(j)] = entity;
-                }
-            }
-            else
-            {
-                vertexToEntity[(*it)->getVertex(j)] = entity;
-            }
-        }
-    }
-}
-
-void SEQ::freeModels(std::vector<GModel*> &models, GModel *global)
-{
-    for(unsigned int i = 0; i < models.size(); i++)
-    {
-        std::vector<GEntity*> entities;
-        models[i]->getEntities(entities);
-        
-        for(unsigned int j = 0; j < entities.size(); j++)
-        {            
-            switch(entities[j]->dim())
-            {
-                case 3:
-                    models[i]->remove(static_cast<GRegion*>(entities[j]));
-                    break;
-                case 2:
-                    models[i]->remove(static_cast<GFace*>(entities[j]));
-                    break;
-                case 1:
-                    models[i]->remove(static_cast<GEdge*>(entities[j]));
-                    break;
-                case 0:
-                    models[i]->remove(static_cast<GVertex*>(entities[j]));
-                    break;
-                default:
-                    break;
-            }
-        }
-    }
-    
-    if(global == nullptr) return;
-        
-    std::vector<GEntity*> entities;
-    global->getEntities(entities);
-    
-    for(unsigned int j = 0; j < entities.size(); j++)
-    {
-        switch(entities[j]->dim())
-        {
-            case 3:
-                global->remove(static_cast<GRegion*>(entities[j]));
-                break;
-            case 2:
-                global->remove(static_cast<GFace*>(entities[j]));
-                break;
-            case 1:
-                global->remove(static_cast<GEdge*>(entities[j]));
-                break;
-            case 0:
-                global->remove(static_cast<GVertex*>(entities[j]));
-                break;
-            default:
-                break;
-        }
-    }
-}
-
-void SEQ::addPhysical(GModel *newModel, GEntity *newEntity, GModel *oldModel, GEntity *oldEntity, GModel *globalModel, GEntity *globalEntity, 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);
-        
-        globalModel->setPhysicalName(name, globalEntity->dim(), number);
-        globalEntity->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]);
-            
-            globalModel->setPhysicalName(name, globalEntity->dim(), oldPhysical[i]);
-            globalEntity->addPhysicalEntity(oldPhysical[i]);
-        }
-    }
-}
-
-
-
-
-
-
-
diff --git a/Partitioning/topology.h b/Partitioning/topology.h
deleted file mode 100644
index 6bad8f7e40445c8681cad0e8d5aa6d831dd5daa6..0000000000000000000000000000000000000000
--- a/Partitioning/topology.h
+++ /dev/null
@@ -1,41 +0,0 @@
-#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
diff --git a/Partitioning/topology_mpi.cpp b/Partitioning/topology_mpi.cpp
deleted file mode 100644
index 202896fa3d0f46da586ca64f35d197d271d6455a..0000000000000000000000000000000000000000
--- a/Partitioning/topology_mpi.cpp
+++ /dev/null
@@ -1,509 +0,0 @@
-#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;
-}
diff --git a/Partitioning/topology_mpi.h b/Partitioning/topology_mpi.h
deleted file mode 100644
index 265573830ab787950b6ba4e53e4baad70876297d..0000000000000000000000000000000000000000
--- a/Partitioning/topology_mpi.h
+++ /dev/null
@@ -1,21 +0,0 @@
-#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