diff --git a/Partitioning/CMakeLists.txt b/Partitioning/CMakeLists.txt
new file mode 100755
index 0000000000000000000000000000000000000000..d69d88ea4334f2040a66cc6dbdd3d53567969e82
--- /dev/null
+++ b/Partitioning/CMakeLists.txt
@@ -0,0 +1,51 @@
+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
new file mode 100755
index 0000000000000000000000000000000000000000..30432292e6d7db5a7688cb0f266b0f01ddfe4a4c
--- /dev/null
+++ b/Partitioning/graph.cpp
@@ -0,0 +1,160 @@
+#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
new file mode 100755
index 0000000000000000000000000000000000000000..d57706f7c225c3b39cbc4a011f2adb2c5aabff2b
--- /dev/null
+++ b/Partitioning/graph.h
@@ -0,0 +1,19 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..7021a881c2ddd525fd332b865512944ae77a17a7
--- /dev/null
+++ b/Partitioning/io.cpp
@@ -0,0 +1,374 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..1dfa6c9319e055cf650a3e550abc917f2183070d
--- /dev/null
+++ b/Partitioning/io.h
@@ -0,0 +1,11 @@
+#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
new file mode 100755
index 0000000000000000000000000000000000000000..3c943c997c7f3f0438652866f3e0a9a86a03a6b2
--- /dev/null
+++ b/Partitioning/main.cpp
@@ -0,0 +1,233 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..c3b75dd73e30e48b18e4504c530f732623a0105e
--- /dev/null
+++ b/Partitioning/main_mpi.cpp
@@ -0,0 +1,200 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..15da23ac2dffb59c740e7022da527d211b8ed88f
--- /dev/null
+++ b/Partitioning/main_mpi.h
@@ -0,0 +1,8 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..656d42f93921ba5ca5e28c1cdfa99c5fcf0fd48c
--- /dev/null
+++ b/Partitioning/topology.cpp
@@ -0,0 +1,1190 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..6bad8f7e40445c8681cad0e8d5aa6d831dd5daa6
--- /dev/null
+++ b/Partitioning/topology.h
@@ -0,0 +1,41 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..202896fa3d0f46da586ca64f35d197d271d6455a
--- /dev/null
+++ b/Partitioning/topology_mpi.cpp
@@ -0,0 +1,509 @@
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..265573830ab787950b6ba4e53e4baad70876297d
--- /dev/null
+++ b/Partitioning/topology_mpi.h
@@ -0,0 +1,21 @@
+#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