Select Git revision
meshPartition.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Anthony Royer authoredAnthony Royer authored
meshPartition.cpp 90.62 KiB
// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.
//
// Contributed by Anthony Royer.
#include <vector>
#include <set>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <limits>
#include <stack>
#include <cstdlib>
#include <map>
#include <utility>
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "GModel.h"
// TODO C++11 remove this nasty stuff
#if __cplusplus >= 201103L
#include <unordered_map>
#define hashmap std::unordered_map
#define hashmapface std::unordered_map\
<MFace, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
Hash_Face, Equal_Face>
#define hashmapedge std::unordered_map\
<MEdge, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
Hash_Edge, Equal_Edge>
#define hashmapvertex std::unordered_map\
<MVertex*, std::vector<std::pair<MElement*, std::vector<unsigned int> > > >
#else
#define hashmap std::map
#define hashmapface std::map\
<MFace, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
Less_Face>
#define hashmapedge std::map\
<MEdge, std::vector<std::pair<MElement*, std::vector<unsigned int> > >,\
Less_Edge>
#define hashmapvertex std::map\
<MVertex*, std::vector<std::pair<MElement*, std::vector<unsigned int> > > >
#endif
#if defined(HAVE_METIS)
#include "OS.h"
#include "Context.h"
#include "partitionRegion.h"
#include "partitionFace.h"
#include "partitionEdge.h"
#include "partitionVertex.h"
#include "ghostRegion.h"
#include "ghostFace.h"
#include "ghostEdge.h"
#include "MFaceHash.h"
#include "MEdgeHash.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "MTrihedron.h"
#include "MElementCut.h"
#include "MPoint.h"
extern "C" {
#include <metis.h>}
// Graph of the mesh for partitioning purposes.
class Graph
{
private:
// The GModel
GModel *_model;
// The number of partitions
unsigned int _nparts;
// The number of elements
size_t _ne;
// The number of nodes
size_t _nn;
// The dimension of the mesh
unsigned int _dim;
// The list of nodes belonging to the ith element of the mesh are stored in
// consecutive locations of eind starting at position eptr[i] up to (but not
// including) position eptr[i+1]. The size of the eind array is of size equal
// to the sum of the number of nodes in all the elements of the mesh.
std::vector<size_t> _eind;
// The size of the eptr array is n + 1, where n is the number of elements in
// the mesh.
std::vector<size_t> _eptr;
// The metis graph structure
std::vector<size_t> _xadj, _adjncy;
// Elements corresponding to each graph elements in eptr
std::vector<MElement*> _element;
// Vertices corresponding to each graph vertices in eptr
std::vector<int> _vertex;
// The width associated to each elements of eptr
std::vector<unsigned int> _vwgt;
// The partitions output from the partitioner
std::vector<unsigned int> _partition;
public:
Graph(GModel * const model)
: _model(model), _nparts(0), _ne(0), _nn(0), _dim(0), _eind(), _eptr(),
_xadj(), _adjncy(), _element(), _vertex(), _vwgt(), _partition()
{
}
void fillDefaultWeights()
{
if(CTX::instance()->mesh.partitionTriWeight == 1 &&
CTX::instance()->mesh.partitionQuaWeight == 1 &&
CTX::instance()->mesh.partitionTetWeight == 1 &&
CTX::instance()->mesh.partitionPyrWeight == 1 &&
CTX::instance()->mesh.partitionPriWeight == 1 &&
CTX::instance()->mesh.partitionHexWeight == 1) return;
_vwgt.resize(_ne);
for(unsigned int i = 0; i < _ne; i++){
if(!_element[i]){
_vwgt[i] = 1;
continue;
}
switch (_element[i]->getType()){
case TYPE_TRI: _vwgt[i] = CTX::instance()->mesh.partitionTriWeight; break;
case TYPE_QUA: _vwgt[i] = CTX::instance()->mesh.partitionQuaWeight; break;
case TYPE_TET: _vwgt[i] = CTX::instance()->mesh.partitionTetWeight; break;
case TYPE_PYR: _vwgt[i] = CTX::instance()->mesh.partitionPyrWeight; break;
case TYPE_PRI: _vwgt[i] = CTX::instance()->mesh.partitionPriWeight; break;
case TYPE_HEX: _vwgt[i] = CTX::instance()->mesh.partitionHexWeight; break;
default: _vwgt[i] = 1; break;
}
}
}
~Graph()
{
clear();
} unsigned int nparts() const { return _nparts; };
size_t ne() const { return _ne; };
size_t nn() const { return _nn; };
unsigned int dim() const { return _dim; };
size_t eind(size_t i) const { return _eind[i]; };
size_t eptr(size_t i) const { return _eptr[i]; };
size_t xadj(size_t i) const { return _xadj[i]; };
std::vector<size_t> &xadj() { return _xadj; };
size_t adjncy(size_t i) const { return _adjncy[i]; };
std::vector<size_t> &adjncy() { return _adjncy; };
MElement* element(size_t i) const { return _element[i]; };
int vertex(size_t i) const { return _vertex[i]; };
std::vector<unsigned int> &vwgt() { return _vwgt; };
unsigned int partition(unsigned int i) const { return _partition[i]; };
std::vector<unsigned int> &partition() { return _partition; };
size_t numNodes() const { return _ne; };
size_t numEdges() const { return _xadj[_ne]/2; };
void nparts(unsigned int nparts) { _nparts = nparts; };
void ne(size_t ne) { _ne = ne; };
void nn(size_t nn) { _nn = nn; };
void dim(unsigned int dim) { _dim = dim; };
void eindResize(size_t size) { _eind.resize(size,0); }
void eind(size_t i, size_t eind) { _eind[i] = eind; };
void eptrResize(size_t size) { _eptr.resize(size, 0); }
void eptr(size_t i, size_t eptr) { _eptr[i] = eptr; };
void elementResize(size_t size){ _element.resize(size, 0); }
void element(size_t i, MElement* element) { _element[i] = element; };
void vertexResize(size_t size) { _vertex.resize(size, -1); }
void adjncy(size_t i, size_t adjncy) { _adjncy[i] = adjncy; };
void vertex(size_t i, int vertex) { _vertex[i] = vertex; };
void vwgt(std::vector<unsigned int> &vwgt) { std::swap(_vwgt, vwgt); };
void partition(std::vector<unsigned int> &partition) { std::swap(_partition, partition); };
void clear()
{
_eind.clear();
_eptr.clear();
_xadj.clear();
_adjncy.clear();
_element.clear();
_vertex.clear();
_vwgt.clear();
_partition.clear();
}
void clearDualGraph()
{
_xadj.clear();
_adjncy.clear();
}
void eraseVertex()
{
for(size_t i = 0; i < _vertex.size(); ++i) _vertex[i] = -1;
}
std::vector< std::set<MElement*> > getBoundaryElements(size_t size = 0)
{
std::vector< std::set<MElement*> > elements
((size ? size : _nparts), std::set<MElement*>());
for(size_t i = 0; i < _ne; ++i){
for(size_t j = _xadj[i]; j < _xadj[i+1]; ++j){
if(_partition[i] != _partition[_adjncy[j]]){
if(_element[i]->getDim() == (int)_dim){
elements[_partition[i]].insert(_element[i]);
}
}
}
}
return elements;
}
void assignGhostCells()
{
std::vector<GEntity*> ghostEntities(_nparts, (GEntity*)NULL);
int elementaryNumber = _model->getMaxElementaryNumber(_dim);
for(unsigned int i = 1; i <= _nparts; ++i){
switch (_dim) {
case 1:
ghostEntities[i - 1] = new ghostEdge(_model, ++elementaryNumber, i);
_model->add(static_cast<ghostEdge*>(ghostEntities[i - 1]));
break;
case 2:
ghostEntities[i - 1] = new ghostFace(_model, ++elementaryNumber, i);
_model->add(static_cast<ghostFace*>(ghostEntities[i - 1]));
break;
case 3:
ghostEntities[i - 1] = new ghostRegion(_model, ++elementaryNumber, i);
_model->add(static_cast<ghostRegion*>(ghostEntities[i - 1]));
break;
default:
break;
}
}
for(size_t i = 0; i < _ne; ++i){
std::set<short> ghostCellsPartition;
for(size_t j = _xadj[i]; j < _xadj[i+1]; j++){
if(_partition[i] != _partition[_adjncy[j]] &&
ghostCellsPartition.find(_partition[_adjncy[j]]) == ghostCellsPartition.end()){
if(_element[i]->getDim() == (int)_dim){
switch (_dim) {
case 1:
static_cast<ghostEdge*>(ghostEntities[_partition[_adjncy[j]]])->
addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
break;
case 2:
static_cast<ghostFace*>(ghostEntities[_partition[_adjncy[j]]])->
addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
break;
case 3:
static_cast<ghostRegion*>(ghostEntities[_partition[_adjncy[j]]])->
addElement(_element[i]->getType(), _element[i], _partition[i] + 1);
break;
default:
break;
}
ghostCellsPartition.insert(_partition[_adjncy[j]]);
}
}
}
}
}
void createDualGraph(bool connectedAll)
{
std::vector<size_t> nptr(_nn+1, 0);
std::vector<size_t> nind(_eptr[_ne], 0);
for(size_t i = 0; i < _ne; ++i){
for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
nptr[_eind[j]]++;
}
}
for(size_t i = 1; i < _nn; ++i) nptr[i] += nptr[i-1];
for(size_t i = _nn; i > 0; --i) nptr[i] = nptr[i-1];
nptr[0] = 0;
for(size_t i = 0; i < _ne; ++i){
for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){ nind[nptr[_eind[j]]++] = i;
}
}
for(unsigned int i = _nn; i > 0; --i) nptr[i] = nptr[i-1];
nptr[0] = 0;
_xadj.resize(_ne+1, 0);
std::vector<size_t> nbrs(_ne, 0);
std::vector<unsigned int> marker(_ne, 0);
for(size_t i = 0; i < _ne; ++i){
size_t l = 0;
for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
for(size_t k = nptr[_eind[j]]; k < nptr[_eind[j]+1]; ++k){
if(nind[k] != i){
if(marker[nind[k]] == 0) nbrs[l++] = nind[k];
marker[nind[k]]++;
}
}
}
unsigned int nbrsNeighbors = 0;
for(size_t j = 0; j < l; j++){
if(marker[nbrs[j]] >=
(connectedAll ? 1 :
_element[i]->numCommonNodesInDualGraph(_element[nbrs[j]])))
nbrsNeighbors++;
marker[nbrs[j]] = 0;
nbrs[j] = 0;
}
_xadj[i] = nbrsNeighbors;
}
for(size_t i = 1; i < _ne; ++i) _xadj[i] = _xadj[i] + _xadj[i-1];
for(size_t i = _ne; i > 0; --i) _xadj[i] = _xadj[i-1];
_xadj[0] = 0;
_adjncy.resize(_xadj[_ne], 0);
for(size_t i = 0; i < _ne; ++i){
size_t l = 0;
for(size_t j = _eptr[i]; j < _eptr[i+1]; ++j){
for(size_t k = nptr[_eind[j]]; k < nptr[_eind[j]+1]; ++k){
if(nind[k] != i){
if (marker[nind[k]] == 0) nbrs[l++] = nind[k];
marker[nind[k]]++;
}
}
}
for(size_t j = 0; j < l; ++j){
if(marker[nbrs[j]] >=
(connectedAll ? 1 :
_element[i]->numCommonNodesInDualGraph(_element[nbrs[j]]))){
_adjncy[_xadj[i]] = nbrs[j];
_xadj[i] = _xadj[i]+1;
}
marker[nbrs[j]] = 0;
nbrs[j] = 0;
}
}
for(size_t i = _ne; i > 0; --i) _xadj[i] = _xadj[i-1];
_xadj[0] = 0;
}
};
template <class ITERATOR>
static void fillElementsToNodesMap(Graph &graph, const GEntity *const entity, size_t &eptrIndex, size_t &eindIndex, size_t &numVertex,
ITERATOR it_beg, ITERATOR it_end)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
const int numVertices = (*it)->getNumPrimaryVertices();
graph.element(eptrIndex++, *it);
graph.eptr(eptrIndex, graph.eptr(eptrIndex-1) + numVertices);
for(int i = 0; i < numVertices; i++){
if(graph.vertex((*it)->getVertex(i)->getNum()-1) == -1){
graph.vertex((*it)->getVertex(i)->getNum()-1, numVertex++);
}
graph.eind(eindIndex, graph.vertex((*it)->getVertex(i)->getNum()-1));
eindIndex++;
}
}
}
static size_t getSizeOfEind(const GModel *const model)
{
size_t size = 0;
// Loop over regions
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
size += 4*(*it)->tetrahedra.size();
size += 8*(*it)->hexahedra.size();
size += 6*(*it)->prisms.size();
size += 5*(*it)->pyramids.size();
size += 4*(*it)->trihedra.size();
}
// Loop over faces
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
size += 3*(*it)->triangles.size();
size += 4*(*it)->quadrangles.size();
}
// Loop over edges
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
size += 2*(*it)->lines.size();
}
// Loop over vertices
for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
size += 1*(*it)->points.size();
}
return size;
}
// Creates a mesh data structure used by Metis routines. Returns: 0 = success, 1
// = no elements found, 2 = error.
static int MakeGraph(GModel *const model, Graph &graph, int selectDim)
{
size_t eindSize = 0;
if(selectDim < 0){
graph.ne(model->getNumMeshElements());
graph.nn(model->getNumMeshVertices());
graph.dim(model->getDim());
graph.elementResize(graph.ne());
graph.vertexResize(model->getMaxVertexNumber());
graph.eptrResize(graph.ne()+1);
graph.eptr(0,0);
eindSize = getSizeOfEind(model);
graph.eindResize(eindSize);
}
else{
GModel* tmp = new GModel();
std::vector<GEntity*> entities;
model->getEntities(entities);
std::set<MVertex*> vertices; for(unsigned int i = 0; i < entities.size(); i++){
if(entities[i]->dim() == selectDim){
switch (entities[i]->dim()) {
case 3:
tmp->add(static_cast<GRegion*>(entities[i]));
break;
case 2:
tmp->add(static_cast<GFace*>(entities[i]));
break;
case 1:
tmp->add(static_cast<GEdge*>(entities[i]));
break;
case 0:
tmp->add(static_cast<GVertex*>(entities[i]));
break;
default:
break;
}
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
for(unsigned int k = 0; k < entities[i]->getMeshElement(j)->getNumVertices() ; k++){
vertices.insert(entities[i]->getMeshElement(j)->getVertex(k));
}
}
}
}
graph.ne(tmp->getNumMeshElements());
graph.nn(vertices.size());
graph.dim(tmp->getDim());
graph.elementResize(graph.ne());
graph.vertexResize(model->getMaxVertexNumber());
graph.eptrResize(graph.ne()+1);
graph.eptr(0,0);
eindSize = getSizeOfEind(tmp);
graph.eindResize(eindSize);
tmp->remove();
delete tmp;
}
size_t eptrIndex = 0;
size_t eindIndex = 0;
size_t numVertex = 0;
if(graph.ne() == 0){
Msg::Error("No mesh elements were found");
return 1;
}
if(graph.dim() == 0){
Msg::Error("Cannot partition a point");
return 1;
}
// Loop over regions
if(selectDim < 0 || selectDim == 3){
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
const GRegion *r = *it;
fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
r->tetrahedra.begin(), r->tetrahedra.end());
fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
r->hexahedra.begin(), r->hexahedra.end());
fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
r->prisms.begin(), r->prisms.end());
fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
r->pyramids.begin(), r->pyramids.end());
fillElementsToNodesMap(graph, r, eptrIndex, eindIndex, numVertex,
r->trihedra.begin(), r->trihedra.end());
}
}
// Loop over faces
if(selectDim < 0 || selectDim == 2){
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
const GFace *f = *it;
fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
f->triangles.begin(), f->triangles.end());
fillElementsToNodesMap(graph, f, eptrIndex, eindIndex, numVertex,
f->quadrangles.begin(), f->quadrangles.end());
}
}
// Loop over edges
if(selectDim < 0 || selectDim == 1){
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
const GEdge *e = *it;
fillElementsToNodesMap(graph, e, eptrIndex, eindIndex, numVertex,
e->lines.begin(), e->lines.end());
}
}
// Loop over vertices
if(selectDim < 0 || selectDim == 0){
for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
GVertex *v = *it;
fillElementsToNodesMap(graph, v, eptrIndex, eindIndex, numVertex,
v->points.begin(), v->points.end());
}
}
return 0;
}
// Partition a graph created by MakeGraph using Metis library. Returns: 0 =
// success, 1 = error, 2 = exception thrown.
static int PartitionGraph(Graph &graph)
{
#ifdef HAVE_METIS
Msg::Info("Running METIS graph partitioner");
try {
int metisOptions[METIS_NOPTIONS];
METIS_SetDefaultOptions((idx_t *)metisOptions);
switch(CTX::instance()->mesh.metisAlgorithm){
case 1: // Recursive
metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_RB;
break;
case 2: // K-way
metisOptions[METIS_OPTION_PTYPE] = METIS_PTYPE_KWAY;
break;
default:
Msg::Info("Unknown partitioning algorithm");
break;
}
switch(CTX::instance()->mesh.metisEdgeMatching){
case 1: // Random matching
metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_RM;
break;
case 2: // Sorted heavy-edge matching
metisOptions[METIS_OPTION_CTYPE] = METIS_CTYPE_SHEM;
break;
default:
Msg::Info("Unknown partitioning edge matching");
break;
}
switch(CTX::instance()->mesh.metisRefinementAlgorithm){ case 1: // FM-based cut refinement
metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_FM;
break;
case 2: // Greedy boundary refinement
metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY;
break;
case 3: // Two-sided node FM refinement
metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP2SIDED;
break;
case 4: // One-sided node FM refinement
metisOptions[METIS_OPTION_RTYPE] = METIS_RTYPE_SEP1SIDED;
break;
default:
Msg::Info("Unknown partitioning refinement algorithm");
break;
}
// C numbering
metisOptions[METIS_OPTION_NUMBERING] = 0;
// Specifies the type of objective
metisOptions[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;
// Forces contiguous partitions.
//metisOptions[METIS_OPTION_CONTIG] = 1;
int objval;
std::vector<unsigned int> epart(graph.ne());
const unsigned int ne = graph.ne();
const int numPart = graph.nparts();
int ncon = 1;
graph.fillDefaultWeights();
int metisError = 0;
graph.createDualGraph(false);
if (metisOptions[METIS_OPTION_PTYPE] == METIS_PTYPE_KWAY){
metisError = METIS_PartGraphKway
((idx_t *)&ne, (idx_t *)&ncon, (idx_t *)graph.xadj().data(),
(idx_t *)graph.adjncy().data(), (idx_t *)graph.vwgt().data(),
(idx_t *)0, 0, (idx_t *)&numPart, 0, 0,
(idx_t *)metisOptions, (idx_t *)&objval, (idx_t *)epart.data());
}
else{
metisError = METIS_PartGraphRecursive
((idx_t *)&ne, (idx_t *)&ncon, (idx_t *)graph.xadj().data(),
(idx_t *)graph.adjncy().data(), (idx_t *)graph.vwgt().data(),
(idx_t *)0, 0, (idx_t *)&numPart, 0, 0,
(idx_t *)metisOptions, (idx_t *)&objval, (idx_t *)epart.data());
}
switch(metisError){
case METIS_OK:
break;
case METIS_ERROR_INPUT:
Msg::Error("METIS input error");
return 1;
case METIS_ERROR_MEMORY:
Msg::Error("METIS memoty error");
return 1;
case METIS_ERROR:
default:
Msg::Error("METIS error");
return 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(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-1){
if(epart[i] != epart[graph.adjncy(j)]){
epart[i] = epart[graph.adjncy(j)];
break;
}
}
}
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-2){
if(epart[i] != epart[graph.adjncy(j)]){
epart[i] = epart[graph.adjncy(j)];
break;
}
}
}
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()-3){
if(epart[i] != epart[graph.adjncy(j)]){
epart[i] = epart[graph.adjncy(j)];
break;
}
}
}
}
graph.partition(epart);
Msg::Info("%d partitions, %d total edge-cuts", numPart, objval);
}
catch(...) {
Msg::Error("METIS exception");
return 2;
}
#endif
return 0;
}
template <class ITERATOR>
static void assignElementsToEntities(GModel *const model,
hashmap<MElement*, unsigned int> &elmToPartition,
std::vector<partitionRegion *> &newRegions,
ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
const unsigned int partition = elmToPartition[(*it)] - 1;
if(!newRegions[partition]){
std::vector<unsigned int> partitions;
partitions.push_back(partition+1);
partitionRegion *dr = new partitionRegion(model, ++elementaryNumber, partitions);
model->add(dr);
newRegions[partition] = dr;
}
newRegions[partition]->addElement((*it)->getType(), *it);
}
}
template <class ITERATOR>
static void assignElementsToEntities(GModel *const model,
hashmap<MElement*, unsigned int> &elmToPartition,
std::vector<partitionFace *> &newFaces,
ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
const unsigned int partition = elmToPartition[(*it)] - 1;
if(!newFaces[partition]){ std::vector<unsigned int> partitions;
partitions.push_back(partition+1);
partitionFace *df = new partitionFace(model, ++elementaryNumber, partitions);
model->add(df);
newFaces[partition] = df;
}
newFaces[partition]->addElement((*it)->getType(), *it);
}
}
template <class ITERATOR>
static void assignElementsToEntities(GModel *const model,
hashmap<MElement*, unsigned int> &elmToPartition,
std::vector<partitionEdge *> &newEdges,
ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
const unsigned int partition = elmToPartition[(*it)] - 1;
if(!newEdges[partition]){
std::vector<unsigned int> partitions;
partitions.push_back(partition+1);
partitionEdge *de = new partitionEdge(model, ++elementaryNumber, 0, 0, partitions);
model->add(de);
newEdges[partition] = de;
}
newEdges[partition]->addElement((*it)->getType(), *it);
}
}
template <class ITERATOR>
static void assignElementsToEntities(GModel *const model,
hashmap<MElement*, unsigned int> &elmToPartition,
std::vector<partitionVertex *> &newVertices,
ITERATOR it_beg, ITERATOR it_end, int &elementaryNumber)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
const unsigned int partition = elmToPartition[(*it)] - 1;
if(!newVertices[partition]){
std::vector<unsigned int> partitions;
partitions.push_back(partition+1);
partitionVertex *dv = new partitionVertex(model, ++elementaryNumber, partitions);
model->add(dv);
newVertices[partition] = dv;
}
newVertices[partition]->addElement((*it)->getType(), *it);
}
}
template <class ITERATOR>
void setVerticesToEntity(GEntity *const entity, ITERATOR it_beg, ITERATOR it_end)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
if(!(*it)->getVertex(i)->onWhat()){
(*it)->getVertex(i)->setEntity(entity);
entity->addMeshVertex((*it)->getVertex(i));
}
}
}
}
// Assign the vertices to its corresponding entity
static void AssignMeshVertices(GModel *model)
{
// Loop over vertices for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
(*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
}
}
(*it)->mesh_vertices.clear();
}
// Loop over edges
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
(*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
}
}
(*it)->mesh_vertices.clear();
}
// Loop over faces
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
(*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
}
}
(*it)->mesh_vertices.clear();
}
// Loop over regions
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++){
for(GModel::size_type j = 0; j < (*it)->getMeshElement(i)->getNumVertices(); j++){
(*it)->getMeshElement(i)->getVertex(j)->setEntity(0);
}
}
(*it)->mesh_vertices.clear();
}
// Loop over vertices
for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
setVerticesToEntity(*it, (*it)->points.begin(), (*it)->points.end());
}
// Loop over edges
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
setVerticesToEntity(*it, (*it)->lines.begin(), (*it)->lines.end());
}
// Loop over faces
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
setVerticesToEntity(*it, (*it)->triangles.begin(), (*it)->triangles.end());
setVerticesToEntity(*it, (*it)->quadrangles.begin(), (*it)->quadrangles.end());
}
// Loop over regions
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
setVerticesToEntity(*it, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
setVerticesToEntity(*it, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
setVerticesToEntity(*it, (*it)->prisms.begin(), (*it)->prisms.end());
setVerticesToEntity(*it, (*it)->pyramids.begin(), (*it)->pyramids.end());
setVerticesToEntity(*it, (*it)->trihedra.begin(), (*it)->trihedra.end());
}
}
static void fillConnectedElements(std::vector< std::set<MElement*> > &connectedElements,
Graph &graph)
{
std::stack<unsigned int> elementStack; std::set<MElement*> elements;
unsigned int startElement = 0;
bool stop = true;
unsigned int size = 0;
int isolatedElements = 0;
do {
//Inititalization
elementStack.push(startElement);
elements.insert(graph.element(startElement));
while(elementStack.size() != 0){
unsigned int top = elementStack.top();
elementStack.pop();
elements.insert(graph.element(top));
for(unsigned int i = graph.xadj(top); i < graph.xadj(top+1); i++){
if(graph.adjncy(i) != 0){
elementStack.push(graph.adjncy(i));
graph.adjncy(i,0);
}
}
}
connectedElements.push_back(elements);
size += elements.size();
elements.clear();
stop = (size == graph.ne() ? true : false);
startElement = 0;
if(!stop){
for(unsigned int i = 0; i < graph.ne(); i++){
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.adjncy(j) != 0){
startElement = i;
i = graph.ne();
break;
}
}
}
if(startElement == 0){
int skipIsolatedElements = 0;
for(unsigned int i = 1; i < graph.ne(); i++){
if(graph.xadj(i) == graph.xadj(i+1)){
if(skipIsolatedElements == isolatedElements){
startElement = i;
isolatedElements++;
break;
}
skipIsolatedElements++;
}
}
}
}
//break;
} while(!stop);
}
static bool dividedNonConnectedEntities(GModel *const model, int dim,
std::set<GRegion*, GEntityLessThan> ®ions,
std::set<GFace*, GEntityLessThan> &faces,
std::set<GEdge*, GEntityLessThan> &edges,
std::set<GVertex*, GEntityLessThan> &vertices)
{
bool ret = false;
// Loop over vertices
if(dim < 0 || dim == 0){
int elementaryNumber = model->getMaxElementaryNumber(0);
for(GModel::const_viter it = vertices.begin(); it != vertices.end(); ++it){ if((*it)->geomType() == GEntity::PartitionVertex){
partitionVertex* vertex = static_cast<partitionVertex*>(*it);
if(vertex->getNumMeshElements() > 1){
ret = true;
for(unsigned int i = 0; i < vertex->getNumMeshElements(); i++){
// Create the new partitionVertex
partitionVertex *pvertex = new partitionVertex
(model, ++elementaryNumber, vertex->getPartitions());
// Assign physicals and parent entity
std::vector<int> physicalTags = vertex->getPhysicalEntities();
for(unsigned int j = 0; j < physicalTags.size(); j++){
pvertex->addPhysicalEntity(physicalTags[j]);
}
pvertex->setParentEntity(vertex->getParentEntity());
// Add to model
model->add(pvertex);
// Add elements
pvertex->addElement(vertex->getMeshElement(i)->getType(),
vertex->getMeshElement(i));
// Move B-Rep
std::vector<GEdge*> BRepEdges = vertex->edges();
if(!BRepEdges.empty()){
for(std::vector<GEdge*>::iterator itBRep = BRepEdges.begin();
itBRep != BRepEdges.end(); ++itBRep){
if(vertex == (*itBRep)->getBeginVertex()){
(*itBRep)->setVertex(pvertex, 1);
pvertex->addEdge(*itBRep);
}
if(vertex == (*itBRep)->getEndVertex()){
(*itBRep)->setVertex(pvertex, -1);
pvertex->addEdge(*itBRep);
}
}
}
}
model->remove(vertex);
vertex->points.clear();
vertex->mesh_vertices.clear();
delete vertex;
}
}
}
}
// Loop over edges
if(dim < 0 || dim == 1){
// We build a graph
Graph graph(model);
graph.ne(model->getNumMeshElements(1));
graph.nn(model->getNumMeshVertices(1));
graph.dim(model->getDim());
graph.elementResize(graph.ne());
graph.vertexResize(model->getMaxVertexNumber());
graph.eptrResize(graph.ne()+1);
graph.eptr(0,0);
const int eindSize = getSizeOfEind(model);
graph.eindResize(eindSize);
int elementaryNumber = model->getMaxElementaryNumber(1);
for(GModel::const_eiter it = edges.begin(); it != edges.end(); ++it){
if((*it)->geomType() == GEntity::PartitionCurve){
partitionEdge* edge = static_cast<partitionEdge*>(*it);
graph.ne(edge->getNumMeshElements());
graph.dim(1);
graph.eptr(0,0);
graph.clearDualGraph(); graph.eraseVertex();
size_t eptrIndex = 0;
size_t eindIndex = 0;
size_t numVertex = 0;
fillElementsToNodesMap(graph, edge, eptrIndex, eindIndex, numVertex,
edge->lines.begin(), edge->lines.end());
graph.nn(numVertex);
graph.createDualGraph(false);
// if a graph contains more than ((n-1)*(n-2))/2 edges
// (where n is the number of nodes), then it is connected.
if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
continue;
}
std::vector< std::set<MElement*> > connectedElements;
fillConnectedElements(connectedElements, graph);
if(connectedElements.size() > 1){
ret = true;
std::vector<GFace*> BRepFaces = edge->faces();
std::vector<int> oldOrientations;
oldOrientations.reserve(BRepFaces.size());
if(!BRepFaces.empty()){
for(std::vector<GFace*>::iterator itBRep = BRepFaces.begin();
itBRep != BRepFaces.end(); ++itBRep){
oldOrientations.push_back((*itBRep)->delEdge(edge));
}
}
for(unsigned int i = 0; i < connectedElements.size(); i++){
// Create the new partitionEdge
partitionEdge *pedge = new partitionEdge
(model, ++elementaryNumber, 0, 0, edge->getPartitions());
// Assign physicals and parent entity
std::vector<int> physicalTags = edge->getPhysicalEntities();
for(unsigned int j = 0; j < physicalTags.size(); j++){
pedge->addPhysicalEntity(physicalTags[j]);
}
pedge->setParentEntity(edge->getParentEntity());
// Add to model
model->add(pedge);
for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
itSet != connectedElements[i].end(); ++itSet){
// Add elements
pedge->addElement((*itSet)->getType(), (*itSet));
}
// Move B-Rep
if(BRepFaces.size() > 0){
int i = 0;
for(std::vector<GFace*>::iterator itBRep = BRepFaces.begin();
itBRep != BRepFaces.end(); ++itBRep){
(*itBRep)->setEdge(pedge, oldOrientations[i]);
pedge->addFace(*itBRep);
i++;
}
}
}
model->remove(edge);
edge->lines.clear();
edge->mesh_vertices.clear();
delete edge;
}
connectedElements.clear(); }
}
}
// Loop over faces
if(dim < 0 || dim == 2){
// We build a graph
Graph graph(model);
graph.ne(model->getNumMeshElements(2));
graph.nn(model->getNumMeshVertices(2));
graph.dim(model->getDim());
graph.elementResize(graph.ne());
graph.vertexResize(model->getMaxVertexNumber());
graph.eptrResize(graph.ne()+1);
graph.eptr(0,0);
const int eindSize = getSizeOfEind(model);
graph.eindResize(eindSize);
int elementaryNumber = model->getMaxElementaryNumber(2);
for(GModel::const_fiter it = faces.begin(); it != faces.end(); ++it){
if((*it)->geomType() == GEntity::PartitionSurface){
partitionFace* face = static_cast<partitionFace*>(*it);
graph.ne(face->getNumMeshElements());
graph.dim(2);
graph.eptr(0,0);
graph.clearDualGraph();
graph.eraseVertex();
size_t eptrIndex = 0;
size_t eindIndex = 0;
size_t numVertex = 0;
fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
face->triangles.begin(), face->triangles.end());
fillElementsToNodesMap(graph, face, eptrIndex, eindIndex, numVertex,
face->quadrangles.begin(), face->quadrangles.end());
graph.nn(numVertex);
graph.createDualGraph(false);
// if a graph contains more than ((n-1)*(n-2))/2 edges
// (where n is the number of nodes), then it is connected.
if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
continue;
}
std::vector< std::set<MElement*> > connectedElements;
fillConnectedElements(connectedElements, graph);
if(connectedElements.size() > 1){
ret = true;
std::list<GRegion*> BRepRegions = face->regions();
std::vector<int> oldOrientations;
if(BRepRegions.size() > 0){
for(std::list<GRegion*>::iterator itBRep = BRepRegions.begin();
itBRep != BRepRegions.end(); ++itBRep){
oldOrientations.push_back((*itBRep)->delFace(face));
}
}
for(unsigned int i = 0; i < connectedElements.size(); i++){
// Create the new partitionFace
partitionFace *pface = new partitionFace
(model, ++elementaryNumber, face->getPartitions());
// Assign physicals and parent entity
std::vector<int> physicalTags = face->getPhysicalEntities();
for(unsigned int j = 0; j < physicalTags.size(); j++){
pface->addPhysicalEntity(physicalTags[j]);
} pface->setParentEntity(face->getParentEntity());
// Add to model
model->add(pface);
for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
itSet != connectedElements[i].end(); ++itSet){
// Add elements
pface->addElement((*itSet)->getType(), (*itSet));
}
// Move B-Rep
if(BRepRegions.size() > 0){
int i = 0;
for(std::list<GRegion*>::iterator itBRep = BRepRegions.begin();
itBRep != BRepRegions.end(); ++itBRep){
(*itBRep)->setFace(pface, oldOrientations[i]);
pface->addRegion(*itBRep);
i++;
}
}
}
model->remove(face);
face->triangles.clear();
face->quadrangles.clear();
face->mesh_vertices.clear();
delete face;
}
connectedElements.clear();
}
}
}
// Loop over regions
if(dim < 0 || dim == 3){
// We build a graph
Graph graph(model);
graph.ne(model->getNumMeshElements(3));
graph.nn(model->getNumMeshVertices(3));
graph.dim(model->getDim());
graph.elementResize(graph.ne());
graph.vertexResize(model->getMaxVertexNumber());
graph.eptrResize(graph.ne()+1);
graph.eptr(0,0);
const int eindSize = getSizeOfEind(model);
graph.eindResize(eindSize);
int elementaryNumber = model->getMaxElementaryNumber(3);
for(GModel::const_riter it = regions.begin(); it != regions.end(); ++it){
if((*it)->geomType() == GEntity::PartitionVolume){
partitionRegion* region = static_cast<partitionRegion*>(*it);
graph.ne(region->getNumMeshElements());
graph.dim(3);
graph.eptr(0,0);
graph.clearDualGraph();
graph.eraseVertex();
size_t eptrIndex = 0;
size_t eindIndex = 0;
size_t numVertex = 0;
fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
region->tetrahedra.begin(), region->tetrahedra.end());
fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
region->hexahedra.begin(), region->hexahedra.end());
fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
region->prisms.begin(), region->prisms.end());
fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
region->pyramids.begin(), region->pyramids.end()); fillElementsToNodesMap(graph, region, eptrIndex, eindIndex, numVertex,
region->trihedra.begin(), region->trihedra.end());
graph.nn(numVertex);
graph.createDualGraph(false);
// if a graph contains more than ((n-1)*(n-2))/2 edges
// (where n is the number of nodes), then it is connected.
if(((graph.numNodes()-1)*(graph.numNodes()-2))/2 < graph.numEdges()){
continue;
}
std::vector< std::set<MElement*> > connectedElements;
fillConnectedElements(connectedElements, graph);
if(connectedElements.size() > 1){
ret = true;
for(unsigned int i = 0; i < connectedElements.size(); i++){
// Create the new partitionRegion
partitionRegion *pregion = new partitionRegion
(model, ++elementaryNumber, region->getPartitions());
// Assign physicals and parent entity
std::vector<int> physicalTags = (*it)->getPhysicalEntities();
for(unsigned int j = 0; j < physicalTags.size(); j++){
pregion->addPhysicalEntity(physicalTags[j]);
}
pregion->setParentEntity(region->getParentEntity());
// Add to model
model->add(pregion);
for(std::set<MElement*>::iterator itSet = connectedElements[i].begin();
itSet != connectedElements[i].end(); ++itSet){
//Add elements
pregion->addElement((*itSet)->getType(), (*itSet));
}
}
model->remove(region);
region->tetrahedra.clear();
region->hexahedra.clear();
region->prisms.clear();
region->pyramids.clear();
region->trihedra.clear();
region->mesh_vertices.clear();
delete region;
}
connectedElements.clear();
}
}
}
return ret;
}
// Create the new volume entities (omega)
static void CreateNewEntities(GModel *const model,
hashmap<MElement*, unsigned int> &elmToPartition)
{
std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
std::set<GFace*, GEntityLessThan> faces = model->getFaces();
std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
int elementaryNumber = model->getMaxElementaryNumber(0);
for(GModel::const_viter it = vertices.begin(); it != vertices.end(); ++it){
std::vector<partitionVertex *> newVertices(model->getNumPartitions(), 0);
assignElementsToEntities(model, elmToPartition, newVertices,
(*it)->points.begin(), (*it)->points.end(), elementaryNumber);
for(unsigned int i = 0; i < model->getNumPartitions(); i++){ if(newVertices[i]){
static_cast<partitionVertex*>(newVertices[i])->setParentEntity((*it));
std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
for(unsigned int j = 0; j < physicalEntities.size(); j++){
newVertices[i]->addPhysicalEntity(physicalEntities[j]);
}
}
}
(*it)->mesh_vertices.clear();
(*it)->points.clear();
}
elementaryNumber = model->getMaxElementaryNumber(1);
for(GModel::const_eiter it = edges.begin(); it != edges.end(); ++it){
std::vector<partitionEdge *> newEdges(model->getNumPartitions(), 0);
assignElementsToEntities(model, elmToPartition, newEdges,
(*it)->lines.begin(), (*it)->lines.end(), elementaryNumber);
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
if(newEdges[i]){
static_cast<partitionEdge*>(newEdges[i])->setParentEntity(*it);
std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
for(unsigned int j = 0; j < physicalEntities.size(); j++){
newEdges[i]->addPhysicalEntity(physicalEntities[j]);
}
}
}
(*it)->mesh_vertices.clear();
(*it)->lines.clear();
}
elementaryNumber = model->getMaxElementaryNumber(2);
for(GModel::const_fiter it = faces.begin(); it != faces.end(); ++it){
std::vector<partitionFace *> newFaces(model->getNumPartitions(), 0);
assignElementsToEntities
(model, elmToPartition, newFaces,
(*it)->triangles.begin(), (*it)->triangles.end(), elementaryNumber);
assignElementsToEntities
(model, elmToPartition, newFaces,
(*it)->quadrangles.begin(), (*it)->quadrangles.end(), elementaryNumber);
std::list<GRegion*> BRepRegions = (*it)->regions();
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
if(newFaces[i]){
static_cast<partitionFace*>(newFaces[i])->setParentEntity(*it);
std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
for(unsigned int j = 0; j < physicalEntities.size(); j++){
newFaces[i]->addPhysicalEntity(physicalEntities[j]);
}
}
}
(*it)->mesh_vertices.clear();
(*it)->triangles.clear();
(*it)->quadrangles.clear();
}
elementaryNumber = model->getMaxElementaryNumber(3);
for(GModel::const_riter it = regions.begin(); it != regions.end(); ++it){
std::vector<partitionRegion *> newRegions(model->getNumPartitions(), 0);
assignElementsToEntities
(model, elmToPartition, newRegions, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end(), elementaryNumber);
assignElementsToEntities
(model, elmToPartition, newRegions,
(*it)->hexahedra.begin(), (*it)->hexahedra.end(), elementaryNumber);
assignElementsToEntities
(model, elmToPartition, newRegions,
(*it)->prisms.begin(), (*it)->prisms.end(), elementaryNumber);
assignElementsToEntities
(model, elmToPartition, newRegions,
(*it)->pyramids.begin(), (*it)->pyramids.end(), elementaryNumber);
assignElementsToEntities
(model, elmToPartition, newRegions,
(*it)->trihedra.begin(), (*it)->trihedra.end(), elementaryNumber);
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
if(newRegions[i]){
static_cast<partitionRegion*>(newRegions[i])->setParentEntity(*it);
std::vector<int> physicalEntities = (*it)->getPhysicalEntities();
for(unsigned int j = 0; j < physicalEntities.size(); j++){
newRegions[i]->addPhysicalEntity(physicalEntities[j]);
}
}
}
(*it)->mesh_vertices.clear();
(*it)->tetrahedra.clear();
(*it)->hexahedra.clear();
(*it)->prisms.clear();
(*it)->pyramids.clear();
(*it)->trihedra.clear();
}
// If we don't create the partition topology let's just assume that the user
// does not care about multi-connected partitions or partition boundaries.
if(!CTX::instance()->mesh.partitionCreateTopology) return;
regions = model->getRegions();
faces = model->getFaces();
edges = model->getEdges();
vertices = model->getVertices();
dividedNonConnectedEntities(model, -1, regions, faces, edges, vertices);
}
static void fillElementToEntity(GModel *const model,
hashmap<MElement*, GEntity*> &elmToEntity, int dim)
{
// Loop over regions
if(dim < 0 || dim == 3){
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
for(std::vector<MTetrahedron*>::iterator itElm = (*it)->tetrahedra.begin();
itElm != (*it)->tetrahedra.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
for(std::vector<MHexahedron*>::iterator itElm = (*it)->hexahedra.begin();
itElm != (*it)->hexahedra.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
for(std::vector<MPrism*>::iterator itElm = (*it)->prisms.begin();
itElm != (*it)->prisms.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
for(std::vector<MPyramid*>::iterator itElm = (*it)->pyramids.begin();
itElm != (*it)->pyramids.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
for(std::vector<MTrihedron*>::iterator itElm = (*it)->trihedra.begin();
itElm != (*it)->trihedra.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
}
}
// Loop over faces
if(dim < 0 || dim == 2){
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){ for(std::vector<MTriangle*>::iterator itElm = (*it)->triangles.begin();
itElm != (*it)->triangles.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
for(std::vector<MQuadrangle*>::iterator itElm = (*it)->quadrangles.begin();
itElm != (*it)->quadrangles.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
}
}
// Loop over edges
if(dim < 0 || dim == 1){
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
for(std::vector<MLine*>::iterator itElm = (*it)->lines.begin();
itElm != (*it)->lines.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
}
}
// Loop over vertices
if(dim < 0 || dim == 0){
for(GModel::const_viter it = model->firstVertex(); it != model->lastVertex(); ++it){
for(std::vector<MPoint*>::iterator itElm = (*it)->points.begin();
itElm != (*it)->points.end(); ++itElm)
elmToEntity.insert(std::pair<MElement*, GEntity*>(*itElm, *it));
}
}
}
static MElement* getReferenceElement(const std::vector< std::pair<MElement*,
std::vector<unsigned int> > > &elementPairs)
{
unsigned int min = std::numeric_limits<unsigned int>::max();
std::vector< std::pair<MElement*, std::vector<unsigned int> > > minSizeElementPairs;
std::vector< std::pair<MElement*, std::vector<unsigned int> > > minSizeElementPairsTmp;
// Take only the elements having the less partition in commun. For exemple we
// take (1,2) and (3,8) but not (2,5,9) or (1,4,5,7)
for(unsigned int i = 0; i < elementPairs.size(); i++){
if(min > elementPairs[i].second.size()){
minSizeElementPairs.clear();
min = elementPairs[i].second.size();
minSizeElementPairs.push_back(elementPairs[i]);
}
else if(min == elementPairs[i].second.size()){
minSizeElementPairs.push_back(elementPairs[i]);
}
}
// Check if the element separated different partitions
if(minSizeElementPairs.size() == elementPairs.size()){
bool isEqual = true;
for(unsigned int i = 1; i < minSizeElementPairs.size(); i++){
if(minSizeElementPairs[i].second != minSizeElementPairs[0].second){
isEqual = false;
break;
}
}
if(isEqual) return 0;
}
while(minSizeElementPairs.size() > 1){
min = std::numeric_limits<unsigned int>::max();
for(unsigned int i = 0; i < minSizeElementPairs.size(); i++){
// The partition vector is sorted thus we can check only the first element
if(minSizeElementPairs[i].second.size() == 0) return minSizeElementPairs[0].first;
if(min > minSizeElementPairs[i].second[0]){
min = minSizeElementPairs[i].second[0];
}
}
for(unsigned int i = 0; i < minSizeElementPairs.size(); i++){
if(min == minSizeElementPairs[i].second[0]){
minSizeElementPairs[i].second.erase(minSizeElementPairs[i].second.begin());
minSizeElementPairsTmp.push_back(minSizeElementPairs[i]);
}
}
minSizeElementPairs.clear();
#if __cplusplus >= 201103L
minSizeElementPairs = std::move(minSizeElementPairsTmp);
#else
minSizeElementPairs = minSizeElementPairsTmp;
#endif
minSizeElementPairsTmp.clear();
}
return minSizeElementPairs[0].first;
}
static void getPartitionInVector(std::vector<unsigned int> &partitions,
const std::vector< std::pair<MElement*,
std::vector<unsigned int> > > &boundaryPair)
{
for(unsigned int i = 0; i < boundaryPair.size(); i++){
for(unsigned int j = 0; j < boundaryPair[i].second.size(); j++){
if(std::find(partitions.begin(), partitions.end(), boundaryPair[i].second[j]) ==
partitions.end()){
partitions.push_back(boundaryPair[i].second[j]);
}
}
}
std::sort(partitions.begin(), partitions.end());
}
static partitionFace *assignPartitionBoundary(
GModel *const model, MFace &me, MElement *reference,
const std::vector<unsigned int> &partitions,
std::multimap<partitionFace *, GEntity *, Less_partitionFace> &pfaces,
hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
{
partitionFace *newEntity = 0;
partitionFace pf(model, 1, partitions);
std::pair< std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator,
std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator>
ret = pfaces.equal_range(&pf);
partitionFace *ppf = 0;
// Create the new partition entity for the mesh
if(ret.first == ret.second){
// Create new entity and add them to the model
ppf = new partitionFace(model, ++numEntity, partitions);
pfaces.insert(std::pair<partitionFace*, GEntity*>(ppf, elementToEntity[reference]));
model->add(ppf);
newEntity = ppf;
}
else{
for(std::multimap<partitionFace*, GEntity*, Less_partitionFace>::iterator it = ret.first;
it != ret.second; ++it){
if(elementToEntity[reference] == (*it).second){
ppf = (*it).first;
}
}
if(!ppf){
// Create new entity and add them to the model
ppf = new partitionFace(model, ++numEntity, partitions);
pfaces.insert(std::pair<partitionFace*, GEntity*>(ppf, elementToEntity[reference]));
model->add(ppf);
newEntity = ppf;
}
}
int numFace = 0;
for(int i = 0; i < reference->getNumFaces(); i++){
if(reference->getFace(i) == me){
numFace = i;
break;
}
}
if(me.getNumVertices() == 3){
std::vector<MVertex*> verts;
reference->getFaceVertices(numFace, verts);
if(verts.size() == 3){
MTriangle *element = new MTriangle(verts);
ppf->addTriangle(element);
}
else if(verts.size() == 6){
MTriangle6 *element = new MTriangle6(verts);
ppf->addTriangle(element);
}
else{
MTriangleN *element = new MTriangleN(verts, verts[0]->getPolynomialOrder());
ppf->addTriangle(element);
}
}
else if(me.getNumVertices() == 4){
std::vector<MVertex*> verts;
reference->getFaceVertices(numFace, verts);
if(verts.size() == 4){
MQuadrangle *element = new MQuadrangle(verts);
ppf->addQuadrangle(element);
}
else if(verts.size() == 8){
MQuadrangle8 *element = new MQuadrangle8(verts);
ppf->addQuadrangle(element);
}
else if(verts.size() == 9){
MQuadrangle9 *element = new MQuadrangle9(verts);
ppf->addQuadrangle(element);
}
else{
MQuadrangleN *element = new MQuadrangleN(verts, verts[0]->getPolynomialOrder());
ppf->addQuadrangle(element);
}
}
return newEntity;
}
static partitionEdge *assignPartitionBoundary(
GModel *const model, MEdge &me, MElement *reference,
const std::vector<unsigned int> &partitions,
std::multimap<partitionEdge *, GEntity *, Less_partitionEdge> &pedges,
hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
{
partitionEdge *newEntity = 0;
partitionEdge pe(model, 1, 0, 0, partitions);
std::pair< std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator,
std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator>
ret = pedges.equal_range(&pe);
partitionEdge *ppe = 0;
// Create the new partition entity for the mesh
if(ret.first == ret.second){
// Create new entity and add them to the model
ppe = new partitionEdge(model, ++numEntity, 0, 0, partitions);
pedges.insert(std::pair<partitionEdge*, GEntity*>(ppe, elementToEntity[reference]));
model->add(ppe); newEntity = ppe;
}
else{
for(std::multimap<partitionEdge*, GEntity*, Less_partitionEdge>::iterator it = ret.first;
it != ret.second; ++it){
if(elementToEntity[reference] == (*it).second){
ppe = (*it).first;
}
}
if(!ppe){
// Create new entity and add them to the model
ppe = new partitionEdge(model, ++numEntity, 0, 0, partitions);
pedges.insert(std::pair<partitionEdge*, GEntity*>(ppe, elementToEntity[reference]));
model->add(ppe);
newEntity = ppe;
}
}
int numEdge = 0;
for(int i = 0; i < reference->getNumEdges(); i++){
if(reference->getEdge(i) == me){
numEdge = i;
break;
}
}
if(me.getNumVertices() == 2){
std::vector<MVertex*> verts;
reference->getEdgeVertices(numEdge, verts);
if(verts.size() == 2){
MLine *element = new MLine(verts);
ppe->addLine(element);
}
else if(verts.size() == 3){
MLine3 *element = new MLine3(verts);
ppe->addLine(element);
}
else{
MLineN *element = new MLineN(verts);
ppe->addLine(element);
}
}
return newEntity;
}
static partitionVertex *assignPartitionBoundary(
GModel *const model, MVertex *ve, MElement *reference,
const std::vector<unsigned int> &partitions,
std::multimap<partitionVertex *, GEntity *, Less_partitionVertex> &pvertices,
hashmap<MElement *, GEntity *> &elementToEntity, int &numEntity)
{
partitionVertex *newEntity = 0;
partitionVertex pv(model, 1, partitions);
std::pair< std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator,
std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator >
ret = pvertices.equal_range(&pv);
partitionVertex *ppv = 0;
// Create the new partition entity for the mesh
if(ret.first == ret.second){
ppv = new partitionVertex(model, ++numEntity, partitions);
pvertices.insert(std::pair<partitionVertex*, GEntity*>(ppv, elementToEntity[reference]));
model->add(ppv);
newEntity = ppv;
}
else{
for(std::multimap<partitionVertex*, GEntity*, Less_partitionVertex>::iterator it =
ret.first; it != ret.second; ++it){ if(elementToEntity[reference] == (*it).second){
ppv = (*it).first;
}
}
if(!ppv){
// Create new entity and add them to the model
ppv = new partitionVertex(model, ++numEntity, partitions);
pvertices.insert(std::pair<partitionVertex*, GEntity*>
(ppv, elementToEntity[reference]));
model->add(ppv);
newEntity = ppv;
}
}
ppv->addPoint(new MPoint(ve));
return newEntity;
}
static int computeOrientation(MElement* reference, MElement* element)
{
if(element->getDim() == 2)
{
std::vector<MVertex*> vertices;
element->getVertices(vertices);
MFace face = element->getFace(0);
for(int i = 0; i < reference->getNumFaces(); i++){
if(reference->getFace(i) == face){
std::vector<MVertex*> referenceVertices;
reference->getFaceVertices(i, referenceVertices);
if(referenceVertices == vertices) return 1;
else return -1;
}
}
}
else if(element->getDim() == 1)
{
std::vector<MVertex*> vertices;
element->getVertices(vertices);
MEdge face = element->getEdge(0);
for(int i = 0; i < reference->getNumEdges(); i++){
if(reference->getEdge(i) == face){
std::vector<MVertex*> referenceVertices;
reference->getEdgeVertices(i, referenceVertices);
if(referenceVertices == vertices) return 1;
else return -1;
}
}
}
else if(element->getDim() == 0)
{
std::vector<MVertex*> vertices;
element->getVertices(vertices);
std::vector<MVertex*> referenceVertices;
reference->getVertices(referenceVertices);
if(referenceVertices[0] == vertices[0]) return 1;
if(referenceVertices[1] == vertices[0]) return -1;
}
return 0;
}
static void assignBrep(GModel *const model, std::map<GEntity*, MElement*>
&boundaryEntityAndRefElement, GEntity *e)
{
if(e->dim() == 2){ partitionFace* entity = static_cast<partitionFace*>(e);
for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
it != boundaryEntityAndRefElement.end(); ++it){
static_cast<GRegion*>(it->first)->setFace
(entity, computeOrientation(it->second, entity->getMeshElement(0)));
entity->addRegion(static_cast<GRegion*>(it->first));
}
}
else if(e->dim() == 1){
partitionEdge* entity = static_cast<partitionEdge*>(e);
for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
it != boundaryEntityAndRefElement.end(); ++it){
static_cast<GFace*>(it->first)->setEdge
(entity, computeOrientation(it->second, entity->getMeshElement(0)));
entity->addFace(static_cast<GFace*>(it->first));
}
}
else if(e->dim() == 0){
partitionVertex* entity = static_cast<partitionVertex*>(e);
for(std::map<GEntity*, MElement*>::iterator it = boundaryEntityAndRefElement.begin();
it != boundaryEntityAndRefElement.end(); ++it){
static_cast<GEdge*>(it->first)->setVertex
(entity, computeOrientation(it->second, entity->getMeshElement(0)));
entity->addEdge(static_cast<GEdge*>(it->first));
}
}
}
void assignNewEntityBRep(Graph &graph, hashmap<MElement*, GEntity*> &elementToEntity)
{
std::set< std::pair<GEntity*, GEntity*> > brepWithoutOri;
hashmap<GEntity*, std::set<std::pair<int, GEntity*> > > brep;
for(unsigned int i = 0; i < graph.ne(); i++){
MElement *current = graph.element(i);
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++)
{
if(current->getDim() == graph.element(graph.adjncy(j))->getDim()+1){
GEntity *g1 = elementToEntity[current];
GEntity *g2 = elementToEntity[graph.element(graph.adjncy(j))];
if(brepWithoutOri.find(std::pair<GEntity*, GEntity*>(g1,g2)) == brepWithoutOri.end()){
const int ori = computeOrientation(current,graph.element(graph.adjncy(j)));
brepWithoutOri.insert(std::pair<GEntity*, GEntity*>(g1,g2));
brep[g1].insert(std::pair<int, GEntity*>(ori,g2));
}
}
}
}
for(hashmap<GEntity*, std::set<std::pair<int, GEntity*> > >::iterator
it = brep.begin(); it != brep.end(); ++it){
switch (it->first->dim()) {
case 3:
for(std::set<std::pair<int, GEntity*> >::iterator itSet =
it->second.begin(); itSet != it->second.end(); ++itSet){
static_cast<GRegion*>(it->first)->setFace
(static_cast<GFace*>(itSet->second), itSet->first);
static_cast<GFace*>(itSet->second)->addRegion
(static_cast<GRegion*>(it->first));
}
break;
case 2:
for(std::set<std::pair<int, GEntity*> >::iterator itSet =
it->second.begin(); itSet != it->second.end(); ++itSet){
static_cast<GFace*>(it->first)->setEdge
(static_cast<GEdge*>(itSet->second), itSet->first);
static_cast<GEdge*>(itSet->second)->addFace
(static_cast<GFace*>(it->first)); }
break;
case 1:
for(std::set<std::pair<int, GEntity*> >::iterator itSet =
it->second.begin(); itSet != it->second.end(); ++itSet){
static_cast<GEdge*>(it->first)->setVertex
(static_cast<GVertex*>(itSet->second), itSet->first);
static_cast<GVertex*>(itSet->second)->addEdge
(static_cast<GEdge*>(it->first));
}
break;
default:
break;
}
}
}
// Create the new entities between each partitions (sigma and bndSigma).
static void CreatePartitionTopology(GModel *const model,
const std::vector< std::set<MElement*> >
&boundaryElements, Graph &meshGraph)
{
const int meshDim = model->getDim();
hashmap<MElement*, GEntity*> elementToEntity;
fillElementToEntity(model, elementToEntity, -1);
assignNewEntityBRep(meshGraph, elementToEntity);
std::multimap<partitionFace*, GEntity*, Less_partitionFace> pfaces;
std::multimap<partitionEdge*, GEntity*, Less_partitionEdge> pedges;
std::multimap<partitionVertex*, GEntity*, Less_partitionVertex> pvertices;
hashmapface faceToElement;
hashmapedge edgeToElement;
hashmapvertex vertexToElement;
std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
std::set<GFace*, GEntityLessThan> faces = model->getFaces();
std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
if (meshDim >= 3){ // Create partition faces
Msg::Info(" - Creating partition faces");
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
it != boundaryElements[i].end(); ++it){
for(int j = 0; j < (*it)->getNumFaces(); j++){
faceToElement[(*it)->getFace(j)].push_back
(std::pair<MElement*, std::vector<unsigned int> >
(*it, std::vector<unsigned int>(1, i+1)));
}
}
}
int numFaceEntity = model->getMaxElementaryNumber(2);
for(hashmapface::const_iterator it = faceToElement.begin();
it != faceToElement.end(); ++it){
MFace f = it->first;
std::vector<unsigned int> partitions;
getPartitionInVector(partitions, it->second);
if(partitions.size() < 2) continue;
MElement* reference = getReferenceElement(it->second);
if(!reference) continue;
partitionFace *pf = assignPartitionBoundary
(model, f, reference, partitions, pfaces, elementToEntity, numFaceEntity);
if(pf){
std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
for(unsigned int i = 0; i < it->second.size(); i++) boundaryEntityAndRefElement.insert
(std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
it->second[i].first));
assignBrep(model, boundaryEntityAndRefElement, pf);
}
}
faceToElement.clear();
faces = model->getFaces();
dividedNonConnectedEntities(model, 2, regions, faces, edges, vertices);
elementToEntity.clear();
fillElementToEntity(model, elementToEntity, 2);
}
if (meshDim >= 2){ // Create partition edges
Msg::Info(" - Creating partition edges");
if (meshDim == 2){
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
it != boundaryElements[i].end(); ++it){
for(int j = 0; j < (*it)->getNumEdges(); j++){
edgeToElement[(*it)->getEdge(j)].push_back
(std::pair<MElement*, std::vector<unsigned int> >
(*it, std::vector<unsigned int>(1, i+1)));
}
}
}
}
else{
Graph subGraph(model);
MakeGraph(model, subGraph, 2);
subGraph.createDualGraph(false);
std::vector<unsigned int>part(subGraph.ne());
int partIndex = 0;
std::map<unsigned int, std::vector<unsigned int> > mapOfPartitions;
unsigned int mapOfPartitionsTag = 0;
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
if((*it)->geomType() == GEntity::PartitionSurface){
std::vector<unsigned int> partitions =
static_cast<partitionFace*>(*it)->getPartitions();
mapOfPartitions.insert(std::pair<unsigned int, std::vector<unsigned int> >
(mapOfPartitionsTag, partitions));
// Must absolutely be in the same order as in the MakeGraph function
for(std::vector<MTriangle*>::iterator itElm = (*it)->triangles.begin();
itElm != (*it)->triangles.end(); ++itElm)
part[partIndex++] = mapOfPartitionsTag;
for(std::vector<MQuadrangle*>::iterator itElm = (*it)->quadrangles.begin();
itElm != (*it)->quadrangles.end(); ++itElm)
part[partIndex++] = mapOfPartitionsTag;
mapOfPartitionsTag++;
}
}
subGraph.partition(part);
std::vector< std::set<MElement*> > subBoundaryElements =
subGraph.getBoundaryElements(mapOfPartitionsTag);
for(unsigned int i = 0; i < mapOfPartitionsTag; i++){
for(std::set<MElement*>::iterator it = subBoundaryElements[i].begin();
it != subBoundaryElements[i].end(); ++it){
for(int j = 0; j < (*it)->getNumEdges(); j++){
edgeToElement[(*it)->getEdge(j)].push_back
(std::pair<MElement*, std::vector<unsigned int> >
(*it, mapOfPartitions[i]));
}
} }
}
int numEdgeEntity = model->getMaxElementaryNumber(1);
for(hashmapedge::const_iterator it = edgeToElement.begin();
it != edgeToElement.end(); ++it){
MEdge e = it->first;
std::vector<unsigned int> partitions;
getPartitionInVector(partitions, it->second);
if(partitions.size() < 2) continue;
MElement* reference = getReferenceElement(it->second);
if(!reference) continue;
partitionEdge* pe = assignPartitionBoundary
(model, e, reference, partitions, pedges, elementToEntity, numEdgeEntity);
if(pe){
std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
for(unsigned int i = 0; i < it->second.size(); i++){
boundaryEntityAndRefElement.insert
(std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
it->second[i].first));
}
assignBrep(model, boundaryEntityAndRefElement, pe);
}
}
edgeToElement.clear();
edges = model->getEdges();
dividedNonConnectedEntities(model, 1, regions, faces, edges, vertices);
elementToEntity.clear();
fillElementToEntity(model, elementToEntity, 1);
}
if (meshDim >= 1){ // Create partition vertices
Msg::Info(" - Creating partition vertices");
if (meshDim == 1){
for(unsigned int i = 0; i < model->getNumPartitions(); i++){
for(std::set<MElement*>::iterator it = boundaryElements[i].begin();
it != boundaryElements[i].end(); ++it){
for(int j = 0; j < (*it)->getNumPrimaryVertices(); j++){
vertexToElement[(*it)->getVertex(j)].push_back
(std::pair<MElement*, std::vector<unsigned int> >
(*it, std::vector<unsigned int>(1,i+1)));
}
}
}
}
else{
Graph subGraph(model);
MakeGraph(model, subGraph, 1);
subGraph.createDualGraph(false);
std::vector<unsigned int> part(subGraph.ne());
int partIndex = 0;
std::map<unsigned int, std::vector<unsigned int> > mapOfPartitions;
unsigned int mapOfPartitionsTag = 0;
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
if((*it)->geomType() == GEntity::PartitionCurve){
std::vector<unsigned int> partitions =
static_cast<partitionEdge*>(*it)->getPartitions();
mapOfPartitions.insert(std::pair<unsigned int, std::vector<unsigned int> >
(mapOfPartitionsTag, partitions));
// Must absolutely be in the same order as in the MakeGraph function
for(std::vector<MLine*>::iterator itElm = (*it)->lines.begin();
itElm != (*it)->lines.end(); ++itElm)
part[partIndex++] = mapOfPartitionsTag;
mapOfPartitionsTag++; }
}
subGraph.partition(part);
std::vector< std::set<MElement*> > subBoundaryElements =
subGraph.getBoundaryElements(mapOfPartitionsTag);
for(unsigned int i = 0; i < mapOfPartitionsTag; i++){
for(std::set<MElement*>::iterator it = subBoundaryElements[i].begin();
it != subBoundaryElements[i].end(); ++it){
for(int j = 0; j < (*it)->getNumPrimaryVertices(); j++){
vertexToElement[(*it)->getVertex(j)].push_back
(std::pair<MElement*, std::vector<unsigned int> >
(*it, mapOfPartitions[i]));
}
}
}
}
int numVertexEntity = model->getMaxElementaryNumber(1);
for(hashmapvertex::const_iterator it = vertexToElement.begin();
it != vertexToElement.end(); ++it){
MVertex *v = it->first;
std::vector<unsigned int> partitions;
getPartitionInVector(partitions, it->second);
if(partitions.size() < 2) continue;
MElement* reference = getReferenceElement(it->second);
if(!reference) continue;
partitionVertex* pv = assignPartitionBoundary
(model, v, reference, partitions, pvertices, elementToEntity, numVertexEntity);
if(pv){
std::map<GEntity*, MElement*> boundaryEntityAndRefElement;
for(unsigned int i = 0; i < it->second.size(); i++)
boundaryEntityAndRefElement.insert
(std::pair<GEntity*, MElement*>(elementToEntity[it->second[i].first],
it->second[i].first));
assignBrep(model, boundaryEntityAndRefElement, pv);
}
}
vertexToElement.clear();
vertices = model->getVertices();
dividedNonConnectedEntities(model, 0, regions, faces, edges, vertices);
}
}
static void addPhysical(GModel *const model, int level,
GEntity *childEntity, hashmap<std::string, int> &nameToNumber,
std::vector<GModel::piter> &iterators, int &numPhysical)
{
if(!CTX::instance()->mesh.partitionCreatePhysicals) return;
unsigned int numPartitions = 0;
if(childEntity->dim() == 3){
numPartitions = static_cast<partitionRegion*>(childEntity)->numPartitions();
}
else if(childEntity->dim() == 2){
numPartitions = static_cast<partitionFace*>(childEntity)->numPartitions();
}
else if(childEntity->dim() == 1){
numPartitions = static_cast<partitionEdge*>(childEntity)->numPartitions();
}
else if(childEntity->dim() == 0){
numPartitions = static_cast<partitionVertex*>(childEntity)->numPartitions();
}
std::string name; if(level == 0){
name = "_omega{";
}
else{
if(numPartitions== 1){
name = "_omicron{";
}
else{
if(level == 1){
name = "_sigma{";
}
else if(level == 2){
name = "_tau{";
}
else if(level == 3){
name = "_upsilon{";
}
}
}
for(unsigned int i = 0; i < numPartitions; i++){
if(i > 0) name += ",";
unsigned int partition = 0;
if(childEntity->dim() == 3){
partition = static_cast<partitionRegion*>(childEntity)->getPartition(i);
}
else if(childEntity->dim() == 2){
partition = static_cast<partitionFace*>(childEntity)->getPartition(i);
}
else if(childEntity->dim() == 1){
partition = static_cast<partitionEdge*>(childEntity)->getPartition(i);
}
else if(childEntity->dim() == 0){
partition = static_cast<partitionVertex*>(childEntity)->getPartition(i);
}
#if __cplusplus >= 201103L
name += std::to_string(partition);
#else
char intToChar[20];
sprintf(intToChar, "%d", partition);
name += intToChar;
#endif
}
name += "}";
int number = 0;
hashmap<std::string, int>::iterator it = nameToNumber.find(name);
if(it == nameToNumber.end()){
number = ++numPhysical;
iterators[childEntity->dim()] = model->setPhysicalName
(iterators[childEntity->dim()], name, childEntity->dim(), number);
nameToNumber.insert(std::pair<std::string, int>(name, number));
}
else{
number = it->second;
}
childEntity->addPhysicalEntity(number);
}
// AssignPhysicalName
static void AssignPhysicalName(GModel *model)
{
std::vector<GModel::piter> iterators;
model->getInnerPhysicalNamesIterators(iterators);
int numPhysical = model->getMaxPhysicalNumber(-1);
hashmap<GEntity*, int> levels;
hashmap<std::string, int> nameToNumber;
// Loop over regions
for(GModel::const_riter it = model->firstRegion(); it != model->lastRegion(); ++it){
if((*it)->geomType() == GEntity::PartitionVolume){ addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(*it, 0));
std::vector<GFace*> listFace = (*it)->faces();
for(std::vector<GFace*>::iterator itF = listFace.begin(); itF != listFace.end(); ++itF){
if(levels.find(*itF) == levels.end()){
const int level = levels[*it]+1;
addPhysical(model, level, *itF, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(*itF, level));
}
}
}
}
// Loop over faces
for(GModel::const_fiter it = model->firstFace(); it != model->lastFace(); ++it){
if((*it)->geomType() == GEntity::PartitionSurface){
if(levels.find(*it) == levels.end()){
addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(*it, 0));
}
std::vector<GEdge*> const& listEdge = (*it)->edges();
for(std::vector<GEdge*>::const_iterator itE = listEdge.begin(); itE != listEdge.end(); ++itE){
if(levels.find(*itE) == levels.end()){
const int level = levels[*it]+1;
addPhysical(model, level, *itE, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(*itE, level));
}
}
}
}
// Loop over edges
for(GModel::const_eiter it = model->firstEdge(); it != model->lastEdge(); ++it){
if((*it)->geomType() == GEntity::PartitionCurve){
if(levels.find(*it) == levels.end()){
addPhysical(model, 0, *it, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(*it, 0));
}
GVertex *v0 = (*it)->getBeginVertex();
if(v0 && levels.find(v0) == levels.end()){
const int level = levels[*it]+1;
addPhysical(model, level, v0, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(v0, level));
}
GVertex *v1 = (*it)->getEndVertex();
if(v1 && levels.find(v1) == levels.end()){
const int level = levels[*it]+1;
addPhysical(model, level, v1, nameToNumber, iterators, numPhysical);
levels.insert(std::pair<GEntity*, int>(v1, level));
}
}
}
}
// Partition a mesh into n parts. Returns: 0 = success, 1 = error
int PartitionMesh(GModel *const model)
{
if(CTX::instance()->mesh.numPartitions <= 0) return 0;
Msg::StatusBar(true, "Partitioning mesh...");
double t1 = Cpu();
Graph graph(model);
if(MakeGraph(model, graph, -1)) return 1;
graph.nparts(CTX::instance()->mesh.numPartitions);
if(PartitionGraph(graph)) return 1;
// Assign partitions to elements
hashmap<MElement*, unsigned int> elmToPartition;
for(unsigned int i = 0; i < graph.ne(); i++){
if(graph.element(i)){
if(graph.nparts() > 1){
elmToPartition.insert(std::pair<MElement*, unsigned int>(graph.element(i),
graph.partition(i)+1));
// Should be removed
graph.element(i)->setPartition(graph.partition(i)+1);
}
else{
elmToPartition.insert(std::pair<MElement*, unsigned int>(graph.element(i), 1));
// Should be removed
graph.element(i)->setPartition(1);
}
}
}
model->setNumPartitions(graph.nparts());
CreateNewEntities(model, elmToPartition);
elmToPartition.clear();
double t2 = Cpu();
Msg::StatusBar(true, "Done partitioning mesh (%g s)", t2 - t1);
if(CTX::instance()->mesh.partitionCreateTopology){
Msg::StatusBar(true, "Creating partition topology...");
std::vector< std::set<MElement*> > boundaryElements = graph.getBoundaryElements();
CreatePartitionTopology(model, boundaryElements, graph);
boundaryElements.clear();
AssignPhysicalName(model);
double t3 = Cpu();
Msg::StatusBar(true, "Done creating partition topology (%g s)", t3 - t2);
}
AssignMeshVertices(model);
if(CTX::instance()->mesh.partitionCreateGhostCells){
graph.clearDualGraph();
graph.createDualGraph(true);
graph.assignGhostCells();
}
return 0;
}
template <class ITERATOR>
static void assignToParent(std::set<MVertex*> &verts, partitionRegion *region,
ITERATOR it_beg, ITERATOR it_end)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
if(region->getParentEntity()->dim() == 3)
region->getParentEntity()->addElement((*it)->getType(), *it);
(*it)->setPartition(0);
for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
if(verts.find((*it)->getVertex(i)) == verts.end()){
(*it)->getVertex(i)->setEntity(region->getParentEntity());
region->getParentEntity()->addMeshVertex((*it)->getVertex(i));
verts.insert((*it)->getVertex(i));
}
}
}
}
template <class ITERATOR>
static void assignToParent(std::set<MVertex*> &verts, partitionFace *face,
ITERATOR it_beg, ITERATOR it_end)
{ for(ITERATOR it = it_beg; it != it_end; ++it){
if(face->getParentEntity()->dim() == 2)
face->getParentEntity()->addElement((*it)->getType(), *it);
(*it)->setPartition(0);
for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
if(verts.find((*it)->getVertex(i)) == verts.end()){
(*it)->getVertex(i)->setEntity(face->getParentEntity());
face->getParentEntity()->addMeshVertex((*it)->getVertex(i));
verts.insert((*it)->getVertex(i));
}
}
}
}
template <class ITERATOR>
static void assignToParent(std::set<MVertex*> &verts, partitionEdge *edge,
ITERATOR it_beg, ITERATOR it_end)
{
for(ITERATOR it = it_beg; it != it_end; ++it){
if(edge->getParentEntity()->dim() == 1)
edge->getParentEntity()->addElement((*it)->getType(), *it);
(*it)->setPartition(0);
for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
if(verts.find((*it)->getVertex(i)) == verts.end()){
(*it)->getVertex(i)->setEntity(edge->getParentEntity());
edge->getParentEntity()->addMeshVertex((*it)->getVertex(i));
verts.insert((*it)->getVertex(i));
}
}
}
}
template <class ITERATOR>
void assignToParent(std::set<MVertex*> &verts, partitionVertex *vertex,
ITERATOR it_beg, ITERATOR it_end)
{
for(ITERATOR it = it_beg; it != it_end; ++it)
{
if(vertex->getParentEntity()->dim() == 0)
vertex->getParentEntity()->addElement((*it)->getType(), *it);
(*it)->setPartition(0);
for(std::size_t i = 0; i < (*it)->getNumVertices(); i++){
if(verts.find((*it)->getVertex(i)) == verts.end()){
(*it)->getVertex(i)->setEntity(vertex->getParentEntity());
vertex->getParentEntity()->addMeshVertex((*it)->getVertex(i));
verts.insert((*it)->getVertex(i));
}
}
}
}
// Un-partition a mesh and return to the initial mesh geomerty. Returns: 0 =
// success, 1 = error.
int UnpartitionMesh(GModel *const model)
{
// make a copy so we can iterate safely (we will remove some entities)
std::set<GRegion*, GEntityLessThan> regions = model->getRegions();
std::set<GFace*, GEntityLessThan> faces = model->getFaces();
std::set<GEdge*, GEntityLessThan> edges = model->getEdges();
std::set<GVertex*, GEntityLessThan> vertices = model->getVertices();
std::set<MVertex*> verts;
// Loop over vertices
for(GModel::viter it = vertices.begin(); it != vertices.end(); ++it){
GVertex *vertex = *it;
if(vertex->geomType() == GEntity::PartitionVertex){
partitionVertex* pvertex = static_cast<partitionVertex*>(vertex);
if(pvertex->getParentEntity()){
assignToParent(verts, pvertex, pvertex->points.begin(), pvertex->points.end());
}
else{
for(unsigned int j = 0; j < pvertex->points.size(); j++)
delete pvertex->points[j];
}
pvertex->points.clear();
pvertex->mesh_vertices.clear();
model->remove(pvertex);
delete pvertex;
}
}
// Loop over edges
for(GModel::eiter it = edges.begin(); it != edges.end(); ++it){
GEdge *edge = *it;
if(edge->geomType() == GEntity::PartitionCurve){
partitionEdge* pedge = static_cast<partitionEdge*>(edge);
if(pedge->getParentEntity()){
assignToParent(verts, pedge, pedge->lines.begin(), pedge->lines.end());
}
else{
for(unsigned int j = 0; j < pedge->lines.size(); j++)
delete pedge->lines[j];
}
pedge->lines.clear();
pedge->mesh_vertices.clear();
pedge->setBeginVertex(0);
pedge->setEndVertex(0);
model->remove(pedge);
delete pedge;
}
else if(edge->geomType() == GEntity::GhostCurve){
model->remove(edge);
delete edge;
}
}
// Loop over faces
for(GModel::fiter it = faces.begin(); it != faces.end(); ++it)
{
GFace *face = *it;
if(face->geomType() == GEntity::PartitionSurface){
partitionFace* pface = static_cast<partitionFace*>(face);
if(pface->getParentEntity()){
assignToParent(verts, pface, pface->triangles.begin(), pface->triangles.end());
assignToParent(verts, pface, pface->quadrangles.begin(), pface->quadrangles.end());
}
else{
for(unsigned int j = 0; j < pface->triangles.size(); j++)
delete pface->triangles[j];
for(unsigned int j = 0; j < pface->quadrangles.size(); j++)
delete pface->quadrangles[j];
}
pface->triangles.clear();
pface->quadrangles.clear();
pface->mesh_vertices.clear();
pface->set(std::vector<GEdge*>());
pface->setOrientations(std::vector<int>());
model->remove(pface);
delete pface;
}
else if(face->geomType() == GEntity::GhostSurface){ model->remove(face);
delete face;
}
}
// Loop over regions
for(GModel::riter it = regions.begin(); it != regions.end(); ++it){
GRegion *region = *it;
if(region->geomType() == GEntity::PartitionVolume){
partitionRegion* pregion = static_cast<partitionRegion*>(region);
if(pregion->getParentEntity()){
assignToParent(verts, pregion, pregion->tetrahedra.begin(), pregion->tetrahedra.end());
assignToParent(verts, pregion, pregion->hexahedra.begin(), pregion->hexahedra.end());
assignToParent(verts, pregion, pregion->prisms.begin(), pregion->prisms.end());
assignToParent(verts, pregion, pregion->pyramids.begin(), pregion->pyramids.end());
assignToParent(verts, pregion, pregion->trihedra.begin(), pregion->trihedra.end());
}
else{
for(unsigned int j = 0; j < pregion->tetrahedra.size(); j++)
delete pregion->tetrahedra[j];
for(unsigned int j = 0; j < pregion->hexahedra.size(); j++)
delete pregion->hexahedra[j];
for(unsigned int j = 0; j < pregion->prisms.size(); j++)
delete pregion->prisms[j];
for(unsigned int j = 0; j < pregion->pyramids.size(); j++)
delete pregion->pyramids[j];
for(unsigned int j = 0; j < pregion->trihedra.size(); j++)
delete pregion->trihedra[j];
}
pregion->tetrahedra.clear();
pregion->hexahedra.clear();
pregion->prisms.clear();
pregion->pyramids.clear();
pregion->trihedra.clear();
pregion->mesh_vertices.clear();
pregion->set(std::vector<GFace*>());
pregion->setOrientations(std::vector<int>());
model->remove(pregion);
delete pregion;
}
else if(region->geomType() == GEntity::GhostVolume){
model->remove(region);
delete region;
}
}
model->setNumPartitions(0);
std::map<std::pair<int, int>, std::string> physicalNames = model->getPhysicalNames();
for(GModel::piter it = physicalNames.begin(); it != physicalNames.end(); ++it){
std::string name = it->second;
if(name[0] == '_'){
model->removePhysicalGroup(it->first.first, it->first.second);
}
}
return 0;
}
// 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){
Graph graph(model);
if(MakeGraph(model, graph, -1)) return 1;
graph.createDualGraph(false);
graph.nparts(npart);
if(elmToPartition.size() != graph.ne()){
Msg::Error("All elements are not partitioned.");
return 1;
}
std::vector<unsigned int> part(graph.ne());
for(unsigned int i = 0; i < graph.ne(); i++){
if(graph.element(i)){
part[i] = elmToPartition[graph.element(i)]-1;
}
}
// Check and correct the topology
// 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(i)->getDim() == graph.element(graph.adjncy(j))->getDim()+1){
if(part[i] != part[graph.adjncy(j)]){
part[i] = part[graph.adjncy(j)];
elmToPartition[graph.element(i)] = part[i]+1;
break;
}
}
}
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->getDim()+2){
if(part[i] != part[graph.adjncy(j)]){
part[i] = part[graph.adjncy(j)];
elmToPartition[graph.element(i)] = part[i]+1;
break;
}
}
}
for(unsigned int j = graph.xadj(i); j < graph.xadj(i+1); j++){
if(graph.element(i)->getDim() == graph.element(graph.adjncy(j))->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);
model->setNumPartitions(graph.nparts());
CreateNewEntities(model, elmToPartition);
elmToPartition.clear();
if(CTX::instance()->mesh.partitionCreateTopology){
Msg::StatusBar(true, "Creating partition topology...");
std::vector< std::set<MElement*> > boundaryElements = graph.getBoundaryElements();
CreatePartitionTopology(model, boundaryElements, graph);
boundaryElements.clear();
AssignPhysicalName(model);
Msg::StatusBar(true, "Done creating partition topology");
}
AssignMeshVertices(model);
if(CTX::instance()->mesh.partitionCreateGhostCells){
graph.clearDualGraph();
graph.createDualGraph(false);
graph.assignGhostCells();
}
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)
{
Msg::Error("Gmsh must be compiled with METIS support to partition meshes");
return 0;
}
int UnpartitionMesh(GModel *const model)
{
return 0;
}
int ConvertOldPartitioningToNewOne(GModel *const model)
{
return 0;
}
int PartitionUsingThisSplit(GModel *const model, unsigned int npart,
hashmap<MElement*, unsigned int> &elmToPartition)
{
return 0;
}
#endif