Commit 842de351 authored by Christophe Geuzaine's avatar Christophe Geuzaine

Merge branch 'simplePartition' into 'master'

Simple partition

See merge request !71
parents 1fa68bcc 157f78df
Pipeline #1093 passed with stage
in 59 minutes and 18 seconds
......@@ -37,6 +37,11 @@ class ghostEdge : public discreteEdge {
bool haveMesh() const { return _haveMesh; }
void haveMesh(bool haveMesh) { _haveMesh = haveMesh; }
virtual std::map<MElement*, unsigned int> &getGhostCells() { return _ghostCells; }
// To make the hidden function visible in ghostEdge
using GEdge::addLine;
using GEdge::addElement;
void addLine(MLine *l, unsigned int onWhichPartition)
{
GEdge::addLine(l);
......
......@@ -41,6 +41,12 @@ class ghostFace : public discreteFace {
bool haveMesh() const { return _haveMesh; }
void haveMesh(bool haveMesh) { _haveMesh = haveMesh; }
virtual std::map<MElement*, unsigned int> &getGhostCells() { return _ghostCells; }
// To make the hidden function visible in ghostFace
using GFace::addTriangle;
using GFace::addQuadrangle;
using GFace::addElement;
void addTriangle(MTriangle *t, unsigned int onWhichPartition)
{
GFace::addTriangle(t);
......
......@@ -45,6 +45,16 @@ class ghostRegion : public discreteRegion {
bool haveMesh() const { return _haveMesh; }
void haveMesh(bool haveMesh) { _haveMesh = haveMesh; }
virtual std::map<MElement*, unsigned int> &getGhostCells() { return _ghostCells; }
// To make the hidden function visible in ghostRegion
using GRegion::addTetrahedron;
using GRegion::addHexahedron;
using GRegion::addPrism;
using GRegion::addPyramid;
using GRegion::addPolyhedron;
using GRegion::addTrihedron;
using GRegion::addElement;
void addTetrahedron(MTetrahedron *t, unsigned int onWhichPartition)
{
GRegion::addTetrahedron(t);
......
......@@ -2822,50 +2822,65 @@ int UnpartitionMesh(GModel *const model)
return 0;
}
// Import a mesh partitionned by a tag given to the element and create the
// topology (omega, sigma, bndSigma, ...). Returns: 0 = success, 1 = no elements
// found.
int ConvertOldPartitioningToNewOne(GModel *const model)
// Create the partition according to the element split given by elmToPartition
// Returns: 0 = success, 1 = no elements found.
int PartitionUsingThisSplit(GModel *const model, unsigned int npart, hashmap<MElement*, unsigned int> &elmToPartition)
{
Msg::StatusBar(true, "Converting old partitioning...");
hashmap<MElement*, unsigned int> elmToPartition;
std::set<unsigned int> partitions;
std::vector<GEntity*> entities;
model->getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(i);
elmToPartition.insert(std::pair<MElement*, unsigned int>(e, e->getPartition()));
partitions.insert(e->getPartition());
}
}
Graph graph(model);
if(MakeGraph(model, graph)) return 1;
createDualGraph(graph, false);
graph.nparts(partitions.size());
graph.nparts(npart);
if(elmToPartition.size() != graph.ne()){
Msg::Error("All elements are not partitioned.");
return 1;
}
unsigned int *part = new unsigned int[graph.ne()];
for(unsigned int i = 0; i < graph.ne(); i++){
if(graph.element(i)){
part[i] = graph.element(i)->getPartition() - 1;
Msg::Info("%d", part[i]);
if(graph.element(i)->getPartition() == 0){
Msg::Error("All elements are not partitioned.");
return 1;
part[i] = elmToPartition[graph.element(i)]-1;
}
}
// Check and correct the topology
for(unsigned int i = 0; i < graph.ne(); i++){
if(graph.element(i)->getDim() == (int)graph.dim()) continue;
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.element(graph.adjncy(j))->getDim() == graph.element(i)->getDim()+1){
if(part[i] != part[graph.adjncy(j)]){
part[i] = part[graph.adjncy(j)];
elmToPartition[graph.element(i)] = part[i]+1;
break;
}
}
if(graph.element(graph.adjncy(j))->getDim() == graph.element(i)->getDim()+2){
if(part[i] != part[graph.adjncy(j)]){
part[i] = part[graph.adjncy(j)];
elmToPartition[graph.element(i)] = part[i]+1;
break;
}
}
if(graph.element(graph.adjncy(j))->getDim() == graph.element(i)->getDim()+3){
if(part[i] != part[graph.adjncy(j)]){
part[i] = part[graph.adjncy(j)];
elmToPartition[graph.element(i)] = part[i]+1;
break;
}
}
}
}
graph.partition(part);
std::vector< std::set<MElement*> > boundaryElements = graph.getBoundaryElements();
model->setNumPartitions(graph.nparts());
CreateNewEntities(model, elmToPartition);
elmToPartition.clear();
Msg::StatusBar(true, "Done converting old partitioning");
if(CTX::instance()->mesh.partitionCreateTopology){
Msg::StatusBar(true, "Creating partition topology...");
CreatePartitionBoundaries(model, boundaryElements);
......@@ -2873,12 +2888,12 @@ int ConvertOldPartitioningToNewOne(GModel *const model)
BuildTopology(model);
Msg::StatusBar(true, "Done creating partition topology");
}
AssignMeshVertices(model);
movePeriodicNodesFromParentToPartitionEntities(model);
if(CTX::instance()->mesh.partitionCreateTopology){
if(CTX::instance()->mesh.partitionCreateGhostCells){
graph.clearDualGraph();
createDualGraph(graph, false);
graph.assignGhostCells();
......@@ -2887,6 +2902,28 @@ int ConvertOldPartitioningToNewOne(GModel *const model)
return 0;
}
// Import a mesh partitionned by a tag given to the element and create the
// topology (omega, sigma, bndSigma, ...). Returns: 0 = success, 1 = no elements
// found.
int ConvertOldPartitioningToNewOne(GModel *const model)
{
Msg::StatusBar(true, "Converting old partitioning...");
hashmap<MElement*, unsigned int> elmToPartition;
std::set<unsigned int> partitions;
std::vector<GEntity*> entities;
model->getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(i);
elmToPartition.insert(std::pair<MElement*, unsigned int>(e, e->getPartition()));
partitions.insert(e->getPartition());
}
}
return PartitionUsingThisSplit(model, partitions.size(), elmToPartition);
}
#else
int PartitionMesh(GModel *const model)
......
......@@ -6,10 +6,18 @@
#ifndef _MESH_PARTITION_H_
#define _MESH_PARTITION_H_
#if __cplusplus >= 201103L
#include <unordered_map>
#define hashmap std::unordered_map
#else
#define hashmap std::map
#endif
class GModel;
int PartitionMesh(GModel *const model);
int UnpartitionMesh(GModel *const model);
int ConvertOldPartitioningToNewOne(GModel *const model);
int PartitionUsingThisSplit(GModel *const model, unsigned int npart, hashmap<MElement*, unsigned int> &elmToPartition);
#endif
......@@ -18,10 +18,18 @@
#include "meshPartition.h"
#endif
#if __cplusplus >= 201103L
#include <unordered_map>
#define hashmap std::unordered_map
#else
#include <map>
#define hashmap std::map
#endif
StringXNumber SimplePartitionOptions_Number[] = {
{GMSH_FULLRC, "NumSlices", NULL, 4.},
{GMSH_FULLRC, "Direction", NULL, 0.},
{GMSH_FULLRC, "CreateBoundaries", NULL, 1.},
{GMSH_FULLRC, "CreateTopology", NULL, 1.},
};
StringXString SimplePartitionOptions_String[] = {
......@@ -41,7 +49,7 @@ std::string GMSH_SimplePartitionPlugin::getHelp() const
return "Plugin(SimplePartition) partitions the current mesh into "
"`NumSlices' slices, along the X-, Y- or Z-axis depending on "
"the value of `Direction' (0,1,2). The plugin creates partition "
"boundaries if `CreateBoundaries' is set.";
"topology if `CreateTopology' is set.";
}
int GMSH_SimplePartitionPlugin::getNbOptions() const
......@@ -64,11 +72,12 @@ StringXString *GMSH_SimplePartitionPlugin::getOptionStr(int iopt)
return &SimplePartitionOptions_String[iopt];
}
PView *GMSH_SimplePartitionPlugin::execute(PView *v)
void GMSH_SimplePartitionPlugin::run()
{
#ifdef HAVE_METIS
int numSlices = (int)SimplePartitionOptions_Number[0].def;
int direction = (int)SimplePartitionOptions_Number[1].def;
int createBoundaries = (int)SimplePartitionOptions_Number[2].def;
int createTopology = (int)SimplePartitionOptions_Number[2].def;
std::vector<std::string> expr(1);
expr[0] = SimplePartitionOptions_String[0].def;
......@@ -84,211 +93,40 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
std::vector<std::string> variables(1);
variables[0] = "t";
mathEvaluator f(expr, variables);
if(expr.empty()) return v;
if(expr.empty()) return;
std::vector<double> values(1), res(1);
for(int p = 0; p <= numSlices; p++){
double t = values[0] = (double)p / (double)numSlices;
if(f.eval(values, res)) t = res[0];
pp[p] = pmin + t * (pmax - pmin);
}
int dim = m->getDim();
std::vector<GEntity*> entities;
m->getEntities(entities);
std::vector<std::set<MFace, Less_Face> > bndFaces(numSlices);
std::vector<std::set<MEdge, Less_Edge> > bndEdges(numSlices);
std::vector<std::set<MElement*> > cutElements(numSlices);
hashmap<MElement*, unsigned int> elmToPartition;
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
if(ge->dim() != dim) continue;
for(unsigned int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
double valmin = pmax;
double valmax = pmin;
bool bnd = false;
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
valmin = std::min(valmin, v->point()[direction]);
valmax = std::max(valmax, v->point()[direction]);
if(v->onWhat() && v->onWhat()->dim() != dim) bnd = true;
}
for(int p = 0; p < numSlices; p++){
if(valmin >= pp[p] && valmin < pp[p + 1]){
e->setPartition(p + 1);
if(bnd){
for(int k = 0; k < e->getNumEdges(); k++)
bndEdges[p].insert(e->getEdge(k));
for(int k = 0; k < e->getNumFaces(); k++)
bndFaces[p].insert(e->getFace(k));
}
if(valmax >= pp[p + 1])
cutElements[p].insert(e);
SPoint3 point = e->barycenter();
for(unsigned int k = 0; k < numSlices; k++){
if(pp[k] < point[direction] && pp[k+1] >= point[direction]){
elmToPartition.insert(std::pair<MElement*, unsigned int>(e, k+1));
break;
}
}
}
}
m->setNumPartitions(numSlices);
// partition lower dimension elements
Msg::Info("Partitioning lower dimension elements");
std::set<MFace, Less_Face> allLowerDimFaces;
std::set<MEdge, Less_Edge> allLowerDimEdges;
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
if(ge->dim() == dim) continue;
for(unsigned int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
if(e->getDim() == 2){
MFace f = e->getFace(0);
if(createBoundaries) allLowerDimFaces.insert(f);
for(int p = 0; p < numSlices; p++){
if(bndFaces[p].find(f) != bndFaces[p].end()){
e->setPartition(p + 1);
break;
}
}
}
else if(e->getDim() == 1){
MEdge f = e->getEdge(0);
if(createBoundaries) allLowerDimEdges.insert(f);
for(int p = 0; p < numSlices; p++){
if(bndEdges[p].find(f) != bndEdges[p].end()){
e->setPartition(p + 1);
break;
}
if(pp[0] == point[direction]){
elmToPartition.insert(std::pair<MElement*, unsigned int>(e, 0+1));
break;
}
}
}
}
opt_mesh_partition_create_topology(0, GMSH_SET, createTopology);
PartitionUsingThisSplit(m, numSlices, elmToPartition);
if(createBoundaries){
Msg::Info("Creating partition boundaries");
#if 0 && defined(HAVE_MESH) // correct & general, but too slow
CreatePartitionBoundaries(m, false, false);
Msg::Info("Renumbering partition boundaries");
std::vector<std::pair<double, GFace*> > parFaces;
std::vector<std::pair<double, GEdge*> > parEdges;
m->getEntities(entities);
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
if(ge->geomType() == GEntity::PartitionSurface ||
ge->geomType() == GEntity::PartitionCurve){
double valmin = pmax;
for(int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
valmin = std::min(valmin, v->point()[direction]);
}
}
if(ge->dim() == 2){
GFace *gf = (GFace*)ge;
parFaces.push_back(std::pair<double, GFace*>(valmin, gf));
m->remove(gf);
}
else{
GEdge *gc = (GEdge*)ge;
parEdges.push_back(std::pair<double, GEdge*>(valmin, gc));
m->remove(gc);
}
}
}
std::sort(parFaces.begin(), parFaces.end());
for(unsigned int i = 0; i < parFaces.size(); i++){
GFace *gf = parFaces[i].second;
gf->setTag(-i-1);
gf->setMeshMaster(-i-1);
m->add(gf);
}
std::sort(parEdges.begin(), parEdges.end());
for(unsigned int i = 0; i < parEdges.size(); i++){
GEdge *ge = parEdges[i].second;
ge->setTag(-i-1);
ge->setMeshMaster(-i-1);
m->add(ge);
}
#else
for(int p = 0; p < numSlices - 1; p++){
std::vector<unsigned int> v2(2);
v2[0] = p + 1;
v2[1] = p + 2;
if(dim == 2){
// create partition edge
partitionEdge *pe = new partitionEdge(m, -(p + 1), 0, 0, v2);
m->add(pe);
// compute boundary of cut surface elements
std::set<MEdge, Less_Edge> myBndEdges;
for(std::set<MElement*>::iterator it = cutElements[p].begin();
it != cutElements[p].end(); it++){
for(int i = 0; i < (*it)->getNumEdges(); i++){
MEdge e = (*it)->getEdge(i);
if(myBndEdges.find(e) == myBndEdges.end())
myBndEdges.insert(e);
else
myBndEdges.erase(e);
}
}
// keep edges whose vertices are all >= than the plane, but are not on a
// curve of the original mesh
for(std::set<MEdge, Less_Edge>::iterator it = myBndEdges.begin();
it != myBndEdges.end(); it++){
bool bnd = true;
for(int j = 0; j < it->getNumVertices(); j++){
if(it->getVertex(j)->point()[direction] < pp[p + 1]){
bnd = false;
break;
}
}
if(bnd){
if(allLowerDimEdges.find(*it) == allLowerDimEdges.end()){
pe->lines.push_back(new MLine(it->getVertex(0), it->getVertex(1)));
}
}
}
}
if(dim == 3){
// create partition face
partitionFace *pf = new partitionFace(m, -(p + 1), v2);
m->add(pf);
// compute boundary of cut elements
std::set<MFace, Less_Face> myBndFaces;
for(std::set<MElement*>::iterator it = cutElements[p].begin();
it != cutElements[p].end(); it++){
for(int i = 0; i < (*it)->getNumFaces(); i++){
MFace f = (*it)->getFace(i);
if(myBndFaces.find(f) == myBndFaces.end())
myBndFaces.insert(f);
else
myBndFaces.erase(f);
}
}
// keep faces whose vertices are all >= than the plane, but are not on a
// surface of the original mesh
for(std::set<MFace, Less_Face>::iterator it = myBndFaces.begin();
it != myBndFaces.end(); it++){
bool bnd = true;
for(int j = 0; j < it->getNumVertices(); j++){
if(it->getVertex(j)->point()[direction] < pp[p + 1]){
bnd = false;
break;
}
}
if(bnd){
if(allLowerDimFaces.find(*it) == allLowerDimFaces.end()){
if (it->getNumVertices() == 3)
pf->triangles.push_back
(new MTriangle(it->getVertex(0), it->getVertex(1), it->getVertex(2)));
else
pf->quadrangles.push_back
(new MQuadrangle(it->getVertex(0), it->getVertex(1),
it->getVertex(2), it->getVertex(3)));
}
}
}
}
}
Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
#endif
}
return v;
}
......@@ -13,7 +13,7 @@ extern "C"
GMSH_Plugin *GMSH_RegisterSimplePartitionPlugin();
}
class GMSH_SimplePartitionPlugin : public GMSH_PostPlugin
class GMSH_SimplePartitionPlugin : public GMSH_MeshPlugin
{
public:
GMSH_SimplePartitionPlugin(){}
......@@ -27,7 +27,7 @@ class GMSH_SimplePartitionPlugin : public GMSH_PostPlugin
StringXNumber* getOption(int iopt);
int getNbOptionsStr() const;
StringXString* getOptionStr(int iopt);
PView *execute(PView *);
void run();
};
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment