Skip to content
Snippets Groups Projects
Commit 930ee39f authored by Anthony Royer's avatar Anthony Royer
Browse files

Delete old files

parent eaff7f64
No related branches found
No related tags found
No related merge requests found
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)
#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
#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
#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;
}
#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 */
#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;
}
#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;
}
#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
This diff is collapsed.
#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
#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;
}
#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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment