Select Git revision
Options.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authored
extended selection buffer to prepare for mesh element picking and the new mesh editing interface
Christophe Geuzaine authoredextended selection buffer to prepare for mesh element picking and the new mesh editing interface
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