Skip to content
Snippets Groups Projects
Commit 735b6ca5 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

ghost cells, take one

parent 8a29a1a7
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GModel.h
+ 6
− 0
View file @ 735b6ca5
......@@ -50,6 +50,10 @@ class GModel
std::vector<MElement*> _elementVectorCache;
std::map<int, MElement*> _elementMapCache;
// ghost cell information (stores partitions for each element acting
// as a ghost cell)
std::multimap<MElement*, short> _ghostCells;
// an octree for fast mesh element lookup
Octree *_octree;
......@@ -314,6 +318,8 @@ class GModel
int getMinPartitionSize() const { return partitionSize[0]; }
int getMaxPartitionSize() const { return partitionSize[1]; }
std::multimap<MElement*, short> &getGhostCells(){ return _ghostCells; }
// perform various coherence tests on the mesh
void checkMeshCoherence(double tolerance);
......
Mesh/meshPartition.cpp
+ 164
− 108
View file @ 735b6ca5
......@@ -75,8 +75,16 @@ template <> struct LFaceTr<MFace>
template <unsigned DIM>
struct MakeGraphFromEntity;
template <unsigned DIM>
struct MatchBoElemToGrVertex;
struct BoElemGr;
typedef std::vector<BoElemGr> BoElemGrVec;
int MakeGraph(GModel *const model, Graph &graph,
BoElemGrVec *const boElemGrVec = 0);
template <unsigned DIM, typename EntIter, typename EntIterBE>
void MakeGraphDIM(const EntIter begin, const EntIter end,
const EntIterBE beginBE, const EntIterBE endBE,
......@@ -104,8 +112,8 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
int PartitionMeshElements( std::vector<MElement*> &elements, meshPartitionOptions &options){
int PartitionMeshElements(std::vector<MElement*> &elements, meshPartitionOptions &options)
{
GModel *tmp_model = new GModel();
GFace *gf = new discreteFace(tmp_model, 1);
std::set<MVertex *> setv;
......@@ -130,7 +138,6 @@ int PartitionMeshElements( std::vector<MElement*> &elements, meshPartitionOption
delete tmp_model;
return 1;
}
int PartitionMeshFace(std::list<GFace*> &cFaces, meshPartitionOptions &options)
......@@ -187,7 +194,7 @@ int PartitionMesh(GModel *const model, meshPartitionOptions &options)
model->recomputeMeshPartitions();
if (options.createPartitionBoundaries)
CreatePartitionBoundaries (model);
CreatePartitionBoundaries (model, true);
Msg::Info("Partitioning complete");
Msg::StatusBar(1, false, "Mesh");
......@@ -733,8 +740,8 @@ void fillit_(std::multimap<MFace,MElement*,Less_Face> &faceToElement,
for(int j=0;j < el->getNumFaces();j++) {
MFace e = el->getFace(j);
faceToElement.insert(std::make_pair(e,el));
}// for every edge of the elem
}// for every elem of the face
}
}
}
template <class ITERATOR>
......@@ -746,8 +753,8 @@ void fillit_(std::multimap<MEdge,MElement*,Less_Edge> &edgeToElement,
for(int j=0;j < el->getNumEdges();j++) {
MEdge e = el->getEdge(j);
edgeToElement.insert(std::make_pair(e,el));
}// for every edge of the elem
}// for every elem of the face
}
}
}
template <class ITERATOR>
......@@ -759,8 +766,8 @@ void fillit_(std::multimap<MVertex*,MElement*> &vertexToElement,
for(int j=0;j < el->getNumVertices();j++) {
MVertex* e = el->getVertex(j);
vertexToElement.insert(std::make_pair(e,el));
}// for every edge of the elem
}// for every elem of the face
}
}
}
void assignPartitionBoundary(GModel *model,
......@@ -948,17 +955,61 @@ void assignPartitionBoundary(GModel *model,
ppv->points.push_back(new MPoint (ve));
}
int CreatePartitionBoundaries(GModel *model)
static void addGhostCells(GEntity *ge,
std::multimap<MVertex*, MElement*> &vertexToElement,
std::multimap<MElement*, short> &ghosts)
{
// get all the nodes on the partition boundary (there are not stored
// so we need to recompute them)
std::set<MVertex*> verts;
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement(i);
for(unsigned int j = 0; j < e->getNumVertices(); j++)
verts.insert(e->getVertex(j));
}
// get all the elements in touching the interface nodes
for(std::set<MVertex*>::iterator it = verts.begin(); it != verts.end(); it++){
MVertex *v = *it;
std::pair<std::multimap<MVertex*, MElement*>::iterator,
std::multimap<MVertex*, MElement*>::iterator> itp =
vertexToElement.equal_range(v);
// get all partitions sharing this node
std::set<short> parts;
for(std::multimap<MVertex*, MElement*>::iterator ite = itp.first;
ite != itp.second; ite++)
parts.insert(ite->second->getPartition());
// update partition info for each element touching the node
for(std::multimap<MVertex*, MElement*>::iterator ite = itp.first;
ite != itp.second; ite++){
MElement *e = ite->second;
for(std::set<short>::iterator its = parts.begin(); its != parts.end(); its++)
if(*its != e->getPartition())
ghosts.insert(std::pair<MElement*, short>(e, *its));
}
}
#if 1
for(std::multimap<MElement*, short>::iterator it = ghosts.begin();
it != ghosts.end(); it++)
printf("ele %d (part %d) ghost %d\n", it->first->getNum(),
it->first->getPartition(), it->second);
#endif
}
int CreatePartitionBoundaries(GModel *model, bool createGhostCells)
{
unsigned numElem[5];
const int meshDim = model->getNumMeshElements(numElem);
std::set<partitionFace*, Less_partitionFace> pfaces;
std::set<partitionEdge*, Less_partitionEdge> pedges;
std::set<partitionVertex*, Less_partitionVertex> pvertices;
std::set<partitionFace*, Less_partitionFace> pfaces;
std::multimap<MFace, MElement*, Less_Face> faceToElement;
std::multimap<MEdge, MElement*, Less_Edge> edgeToElement;
std::multimap<MVertex*, MElement*> vertexToElement;
// assign partition faces
if (meshDim == 3){
std::multimap<MFace,MElement*,Less_Face> faceToElement;
for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it){
fillit_(faceToElement, (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
fillit_(faceToElement, (*it)->hexahedra.begin(), (*it)->hexahedra.end());
......@@ -966,7 +1017,6 @@ int CreatePartitionBoundaries(GModel *model)
fillit_(faceToElement, (*it)->pyramids.begin(), (*it)->pyramids.end());
fillit_(faceToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
{
std::multimap<MFace, MElement*, Less_Face>::iterator it = faceToElement.begin();
Equal_Face oper;
while (it != faceToElement.end()){
......@@ -980,11 +1030,9 @@ int CreatePartitionBoundaries(GModel *model)
assignPartitionBoundary(model, e, pfaces, voe);
}
}
}
// assign partition edges
if (meshDim > 1){
std::multimap<MEdge,MElement*,Less_Edge> edgeToElement;
if (meshDim == 2){
for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
fillit_(edgeToElement, (*it)->triangles.begin(), (*it)->triangles.end());
......@@ -1001,7 +1049,6 @@ int CreatePartitionBoundaries(GModel *model)
fillit_(edgeToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
}
{
std::multimap<MEdge, MElement*, Less_Edge>::iterator it = edgeToElement.begin();
Equal_Edge oper;
while (it != edgeToElement.end()){
......@@ -1016,11 +1063,9 @@ int CreatePartitionBoundaries(GModel *model)
}
//splitBoundaryEdges(model,pedges);
}
}
// make partition vertices
if (meshDim > 1){
std::multimap<MVertex*,MElement*> vertexToElement;
if (meshDim == 2){
for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){
fillit_(vertexToElement, (*it)->triangles.begin(), (*it)->triangles.end());
......@@ -1037,7 +1082,6 @@ int CreatePartitionBoundaries(GModel *model)
fillit_(vertexToElement, (*it)->polyhedra.begin(), (*it)->polyhedra.end());
}
}
{
std::multimap<MVertex*, MElement*>::iterator it = vertexToElement.begin();
while (it != vertexToElement.end()){
MVertex *v = it->first;
......@@ -1050,6 +1094,18 @@ int CreatePartitionBoundaries(GModel *model)
assignPartitionBoundary(model, v, pvertices, voe, pedges, pfaces);
}
}
if(createGhostCells){
std::multimap<MElement*, short> &ghosts(model->getGhostCells());
ghosts.clear();
if(meshDim == 2)
for(std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedges.begin();
it != pedges.end(); it++)
addGhostCells(*it, vertexToElement, ghosts);
else if(meshDim == 3)
for(std::set<partitionFace*, Less_partitionFace>::iterator it = pfaces.begin();
it != pfaces.end(); it++)
addGhostCells(*it, vertexToElement, ghosts);
}
return 1;
......
Mesh/meshPartition.h
+ 5
− 11
View file @ 735b6ca5
......@@ -3,22 +3,18 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _PARTITION_H_
#define _PARTITION_H_
#ifndef _MESH_PARTITION_H_
#define _MESH_PARTITION_H_
#include <vector>
#include "partitionEdge.h"
#include "GFaceCompound.h"
#include "GFace.h"
struct meshPartitionOptions;
struct BoElemGr;
class GModel;
class GFace;
class Graph;
typedef std::vector<BoElemGr> BoElemGrVec;
typedef enum {LAPLACIAN= 0, MULTILEVEL=1} typeOfPartition;
struct meshPartitionOptions;
/*******************************************************************************
*
......@@ -26,14 +22,12 @@ typedef enum {LAPLACIAN= 0, MULTILEVEL=1} typeOfPartition;
*
******************************************************************************/
int MakeGraph(GModel *const model, Graph &graph,
BoElemGrVec *const boElemGrVec = 0);
int PartitionGraph(Graph &graph, meshPartitionOptions &options);
int PartitionMesh(GModel *const model, meshPartitionOptions &options);
int PartitionMeshFace(std::list<GFace*> &cFaces, meshPartitionOptions &options);
int PartitionMeshElements(std::vector<MElement*> &elements, meshPartitionOptions &options);
bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts);
int CreatePartitionBoundaries (GModel *model);
int CreatePartitionBoundaries(GModel *model, bool createGhostCells);
void splitBoundaryEdges(GModel *model, std::set<partitionEdge*, Less_partitionEdge> &newEdges);
void createPartitionFaces(GModel *model, GFaceCompound * gf, int num_parts,
......
Mesh/multiscalePartition.cpp
+ 69
− 79
View file @ 735b6ca5
......@@ -10,8 +10,8 @@
static void recur_connect(MVertex *v,
std::multimap<MVertex*,MEdge> &v2e,
std::set<MEdge,Less_Edge> &group,
std::set<MVertex*> &touched){
std::set<MVertex*> &touched)
{
if (touched.find(v) != touched.end())return;
touched.insert(v);
......@@ -27,7 +27,6 @@ static void recur_connect (MVertex *v,
static int connected_bounds (std::vector<MEdge> &edges, std::vector<std::vector<MEdge> > &boundaries)
{
std::multimap<MVertex*,MEdge> v2e;
for (unsigned i = 0; i < edges.size(); ++i){
for (int j=0;j<edges[i].getNumVertices();j++){
......@@ -50,7 +49,8 @@ static int connected_bounds (std::vector<MEdge> &edges, std::vector<std::vector
static int getGenus (std::vector<MElement *> &elements,
std::vector<std::vector<MEdge> > &boundaries){
std::vector<std::vector<MEdge> > &boundaries)
{
//We suppose MElements are simply connected
......@@ -94,7 +94,6 @@ static int getGenus (std::vector<MElement *> &elements,
static int getAspectRatio(std::vector<MElement *> &elements,
std::vector<std::vector<MEdge> > &boundaries)
{
std::set<MVertex*> vs;
for(unsigned int i = 0; i < elements.size(); i++){
MElement *e = elements[i];
......@@ -116,10 +115,10 @@ static int getAspectRatio (std::vector<MElement *> &elements,
double H = norm(SVector3(bb.max(), bb.min()));
double D = H;
for (unsigned i=0; i< boundaries.size(); i++){
for (unsigned int i = 0; i < boundaries.size(); i++){
std::set<MVertex*> vb;
std::vector<MEdge> iBound = boundaries[i];
for (unsigned j=0; j< iBound.size(); j++){
for (unsigned int j = 0; j < iBound.size(); j++){
MEdge e = iBound[j];
vb.insert(e.getVertex(0));
vb.insert(e.getVertex(1));
......@@ -135,32 +134,27 @@ static int getAspectRatio (std::vector<MElement *> &elements,
int AR = (int)ceil(H/D);
return AR;
}
static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR){
static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR)
{
std::vector<std::vector<MEdge> > boundaries;
genus = getGenus(elements, boundaries);
AR = getAspectRatio(elements, boundaries);
return;
}
static void partitionRegions (std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > &regions){
for (unsigned i=0;i<elements.size();++i){
static void partitionRegions(std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > &regions)
{
for (unsigned int i = 0; i < elements.size(); ++i){
MElement *e = elements[i];
int part = e->getPartition();
regions[part-1].push_back(e);
}
return;
}
static void printLevel ( std::vector<MElement *> &elements, int recur, int region){
static void printLevel(std::vector<MElement *> &elements, int recur, int region)
{
char fn[256];
sprintf(fn, "part_%d_%d.msh", recur, region);
double version = 2.0;
......@@ -178,7 +172,7 @@ static void printLevel ( std::vector<MElement *> &elements, int recur, int regio
fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
fprintf(fp, "$EndMeshFormat\n");
fprintf(fp,"$Nodes\n%d\n",vs.size());
fprintf(fp,"$Nodes\n%d\n", (int)vs.size());
std::set<MVertex*> :: iterator it = vs.begin();
int index = 1;
for (; it != vs.end() ; ++it){
......@@ -188,18 +182,17 @@ static void printLevel ( std::vector<MElement *> &elements, int recur, int regio
}
fprintf(fp,"$EndNodes\n",elements.size());
fprintf(fp,"$Elements\n%d\n",elements.size());
fprintf(fp,"$Elements\n%d\n", (int)elements.size());
for (int i=0;i<elements.size();i++){
elements[i]->writeMSH(fp,version);
}
fprintf(fp,"$EndElements\n%d\n",elements.size());
fprintf(fp,"$EndElements\n%d\n", (int)elements.size());
fclose(fp);
}
multiscalePartition::multiscalePartition(std::vector<MElement *> &elements, int nbParts, typeOfPartition method)
{
options = CTX::instance()->partitionOptions;
options.num_partitions = nbParts;
options.partitioner = 1; //1 CHACO, 2 METIS
......@@ -218,8 +211,8 @@ multiscalePartition::multiscalePartition (std::vector<MElement *> &elements, int
partition(*level, nbParts, method);
totalParts = assembleAllPartitions();
}
void multiscalePartition::setNumberOfPartitions(int &nbParts)
{
options.num_partitions = nbParts;
......@@ -228,8 +221,8 @@ void multiscalePartition:: setNumberOfPartitions(int &nbParts)
}
}
void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfPartition method){
void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfPartition method)
{
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
if (method == LAPLACIAN){
......@@ -284,8 +277,8 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
#endif
}
int multiscalePartition::assembleAllPartitions(){
int multiscalePartition::assembleAllPartitions()
{
int nbParts = 1;
for (unsigned i = 0; i< levels.size(); i++){
......@@ -301,7 +294,4 @@ int multiscalePartition::assembleAllPartitions(){
}
return nbParts - 1;
}
Mesh/multiscalePartition.h
+ 2
− 0
View file @ 735b6ca5
......@@ -21,6 +21,8 @@ struct partitionLevel {
std::vector<MElement *> elements;
};
typedef enum {LAPLACIAN= 0, MULTILEVEL=1} typeOfPartition;
class multiscalePartition{
private:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment