diff --git a/Geo/GEntity.h b/Geo/GEntity.h index 1fe452ccee3945d14d75c3d75ff2a3c5fd88b5a2..323dad686e34ee43877e3254fb552c4af06ad96d 100644 --- a/Geo/GEntity.h +++ b/Geo/GEntity.h @@ -101,7 +101,7 @@ class GEntity { DiscreteVolume, CompoundVolume, PartitionVertex, - PartitionEdge, + PartitionCurve, PartitionSurface }; @@ -144,9 +144,9 @@ class GEntity { "Volume", "Discrete volume", "Compound Volume", - "Partition Vertex", - "Partition Edge", - "Partition Surface" + "Partition vertex", + "Partition curve", + "Partition surface" }; unsigned int type = (unsigned int)geomType(); if(type >= sizeof(name) / sizeof(name[0])) diff --git a/Geo/partitionEdge.h b/Geo/partitionEdge.h new file mode 100644 index 0000000000000000000000000000000000000000..4afe2054033497a3208d5e884646d64f4121dccf --- /dev/null +++ b/Geo/partitionEdge.h @@ -0,0 +1,42 @@ +// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. + +#ifndef _PARTITION_EDGE_H_ +#define _PARTITION_EDGE_H_ + +#include "GModel.h" +#include "GEdge.h" +#include "discreteEdge.h" + +class partitionEdge : public discreteEdge { + public: + std::vector<int> _partitions; + public: + partitionEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1, + std::vector<int> &partitions) + : discreteEdge(model,num,_v0,_v1),_partitions(partitions) + { + std::sort(_partitions.begin(),_partitions.end()); + } + virtual ~partitionEdge() {} + virtual GeomType geomType() const { return PartitionCurve; } +}; + + +struct Less_partitionEdge : + public std::binary_function<partitionEdge*, partitionEdge*, bool> { + bool operator()(const partitionEdge* e1, const partitionEdge* e2) const{ + if (e1->_partitions.size() < e2->_partitions.size())return true; + if (e1->_partitions.size() > e2->_partitions.size())return false; + for (int i=0;i<e1->_partitions.size();i++){ + if (e1->_partitions[i] < e2->_partitions[i])return true; + if (e1->_partitions[i] > e2->_partitions[i])return false; + } + return false; + } +}; + + +#endif diff --git a/Geo/partitionFace.h b/Geo/partitionFace.h new file mode 100644 index 0000000000000000000000000000000000000000..53e860fec0f4ae8a13b38c8d1f91ad9a3ddbca41 --- /dev/null +++ b/Geo/partitionFace.h @@ -0,0 +1,40 @@ +// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. + +#ifndef _PARTITION_FACE_H_ +#define _PARTITION_FACE_H_ + +#include "GModel.h" +#include "discreteFace.h" + +class partitionFace : public discreteFace { + public: + std::vector<int> _partitions; + public: + partitionFace(GModel *model, int num, std::vector<int> &partitions) + : discreteFace(model,num),_partitions(partitions) + { + std::sort(_partitions.begin(),_partitions.end()); + } + virtual ~partitionFace() {} + virtual GeomType geomType() const { return PartitionSurface; } +}; + + +struct Less_partitionFace : + public std::binary_function<partitionFace*, partitionFace*, bool> { + bool operator()(const partitionFace* e1, const partitionFace* e2) const{ + if (e1->_partitions.size() < e2->_partitions.size())return true; + if (e1->_partitions.size() > e2->_partitions.size())return false; + for (int i=0;i<e1->_partitions.size();i++){ + if (e1->_partitions[i] < e2->_partitions[i])return true; + if (e1->_partitions[i] > e2->_partitions[i])return false; + } + return false; + } +}; + + +#endif diff --git a/Geo/partitionVertex.h b/Geo/partitionVertex.h new file mode 100644 index 0000000000000000000000000000000000000000..9427088b6b4189e02a8e54502a102aec37a5e08e --- /dev/null +++ b/Geo/partitionVertex.h @@ -0,0 +1,40 @@ +// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. + +#ifndef _PARTITION_VERTEX_H_ +#define _PARTITION_VERTEX_H_ + +#include "GModel.h" +#include "discreteVertex.h" + +class partitionVertex : public discreteVertex { + public: + std::vector<int> _partitions; + public: + partitionVertex(GModel *model, int num, std::vector<int> &partitions) + : discreteVertex(model,num),_partitions(partitions) + { + std::sort(_partitions.begin(),_partitions.end()); + } + virtual ~partitionVertex() {} + virtual GeomType geomType() const { return PartitionVertex; } +}; + + +struct Less_partitionVertex : + public std::binary_function<partitionVertex*, partitionVertex*, bool> { + bool operator()(const partitionVertex* e1, const partitionVertex* e2) const{ + if (e1->_partitions.size() < e2->_partitions.size())return true; + if (e1->_partitions.size() > e2->_partitions.size())return false; + for (int i=0;i<e1->_partitions.size();i++){ + if (e1->_partitions[i] < e2->_partitions[i])return true; + if (e1->_partitions[i] > e2->_partitions[i])return false; + } + return false; + } +}; + + +#endif diff --git a/Graphics/drawGeom.cpp b/Graphics/drawGeom.cpp index 60d2b979e17e89c7213e8ffee30a00a794c2ef38..cfe11bf5a27d79db6f9173127528a96caf2ccbb3 100644 --- a/Graphics/drawGeom.cpp +++ b/Graphics/drawGeom.cpp @@ -95,6 +95,7 @@ class drawGEdge { { if(!e->getVisibility()) return; if(e->geomType() == GEntity::DiscreteCurve) return; + if(e->geomType() == GEntity::PartitionCurve) return; if(e->geomType() == GEntity::BoundaryLayerCurve) return; bool select = (_ctx->render_mode == drawContext::GMSH_SELECT && diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp index 2aa005e2c812b4c4f548504db78dd814887366a5..e396ef6eec42ba73487ce0a12ba7684f5d8898aa 100644 --- a/Mesh/meshPartition.cpp +++ b/Mesh/meshPartition.cpp @@ -15,8 +15,15 @@ #include "meshPartitionOptions.h" #include "MTriangle.h" #include "MQuadrangle.h" +#include "MTetrahedron.h" +#include "MHexahedron.h" +#include "MPrism.h" +#include "MPyramid.h" #include "MElementCut.h" +#include "MPoint.h" +#include "partitionVertex.h" #include "partitionEdge.h" +#include "partitionFace.h" //--Prototype for Chaco interface @@ -129,7 +136,8 @@ int PartitionMesh(GModel *const model, meshPartitionOptions &options) model->recomputeMeshPartitions(); - if (options.createPartitionBoundaries)CreatePartitionBoundaries (model); + // if (options.createPartitionBoundaries) + CreatePartitionBoundaries (model); Msg::Info("Partitioning complete"); Msg::StatusBar(1, false, "Mesh"); @@ -650,6 +658,19 @@ struct MatchBoElemToGrVertex } }; + +template <class ITERATOR> +void fillit_ ( std::multimap<MFace,MElement*,Less_Face> &faceToElement, + ITERATOR it_beg, ITERATOR it_end){ + for (ITERATOR IT = it_beg; IT != it_end ; ++IT){ + MElement *el = *IT; + 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> void fillit_ ( std::multimap<MEdge,MElement*,Less_Edge> &edgeToElement, ITERATOR it_beg, ITERATOR it_end){ @@ -675,10 +696,48 @@ void fillit_ ( std::multimap<MVertex*,MElement*> &vertexToElement, } +void assignPartitionBoundary (GModel *model, + MFace &me, + std::set<partitionFace*, Less_partitionFace> &pfaces, + std::vector<MElement*> &v){ + + std::vector<int> v2; + v2.push_back(v[0]->getPartition()); + for (int i=1;i<v.size();i++){ + bool found = false; + for (int j=0;j<v2.size();j++){ + if (v[i]->getPartition() == v2[j]){ + found = true; + break; + } + } + if (!found)v2.push_back(v[i]->getPartition()); + } + if (v2.size() < 2)return; + + partitionFace pe (model, 1, v2); + std::set<partitionFace*, Less_partitionFace>::iterator it = pfaces.find(&pe); + partitionFace *ppe; + if (it == pfaces.end()){ + ppe = new partitionFace(model, -pfaces.size()-1, v2); + pfaces.insert (ppe); + model->add(ppe); + printf("created partitionFace %d (",ppe->tag()); + for (int i=0;i<v2.size();++i)printf("%d ",v2[i]); + printf(")\n"); + } + else ppe = *it; + // to finish !! + ppe->triangles.push_back(new MTriangle (me.getVertex(0),me.getVertex(1), + me.getVertex(2))); +} + + void assignPartitionBoundary (GModel *model, MEdge &me, std::set<partitionEdge*, Less_partitionEdge> &pedges, - std::vector<MElement*> &v){ + std::vector<MElement*> &v, + std::set<partitionFace*, Less_partitionFace> &pfaces){ std::vector<int> v2; v2.push_back(v[0]->getPartition()); @@ -693,6 +752,11 @@ void assignPartitionBoundary (GModel *model, if (!found)v2.push_back(v[i]->getPartition()); } if (v2.size() < 2)return; + + + partitionFace pf (model, 1,v2); + std::set<partitionFace*, Less_partitionFace>::iterator itf = pfaces.find(&pf); + if (itf != pfaces.end())return; partitionEdge pe (model, 1, 0, 0, v2); std::set<partitionEdge*, Less_partitionEdge>::iterator it = pedges.find(&pe); @@ -701,35 +765,171 @@ void assignPartitionBoundary (GModel *model, ppe = new partitionEdge(model, -pedges.size()-1, 0, 0, v2); pedges.insert (ppe); model->add(ppe); + printf("created partitionEdge %d (",ppe->tag()); + for (int i=0;i<v2.size();++i)printf("%d ",v2[i]); + printf(")\n"); } else ppe = *it; ppe->lines.push_back(new MLine (me.getVertex(0),me.getVertex(1))); } +void assignPartitionBoundary (GModel *model, + MVertex *ve, + std::set<partitionVertex*, Less_partitionVertex> &pvertices, + std::vector<MElement*> &v, + std::set<partitionEdge*, Less_partitionEdge> &pedges, + std::set<partitionFace*, Less_partitionFace> &pfaces){ + + std::vector<int> v2; + v2.push_back(v[0]->getPartition()); + for (int i=1;i<v.size();i++){ + bool found = false; + for (int j=0;j<v2.size();j++){ + if (v[i]->getPartition() == v2[j]){ + found = true; + break; + } + } + if (!found)v2.push_back(v[i]->getPartition()); + } + if (v2.size() < 2)return; + + partitionFace pf (model, 1,v2); + std::set<partitionFace*, Less_partitionFace>::iterator itf = pfaces.find(&pf); + if (itf != pfaces.end())return; + + partitionEdge pe (model, 1,0,0,v2); + std::set<partitionEdge*, Less_partitionEdge>::iterator ite = pedges.find(&pe); + if (ite != pedges.end())return; + + partitionVertex pv (model, 1,v2); + std::set<partitionVertex*, Less_partitionVertex>::iterator it = pvertices.find(&pv); + partitionVertex *ppv; + if (it == pvertices.end()){ + ppv = new partitionVertex(model, -pvertices.size()-1,v2); + pvertices.insert (ppv); + model->add(ppv); + printf("created partitionVertex %d (",ppv->tag()); + for (int i=0;i<v2.size();++i)printf("%d ",v2[i]); + printf(")\n"); + } + else ppv = *it; + ppv->points.push_back(new MPoint (ve)); +} + + int CreatePartitionBoundaries (GModel *model) { - std::multimap<MEdge,MElement*,Less_Edge> edgeToElement; + + unsigned numElem[5]; + const int meshDim = model->getNumMeshElements(numElem); std::set<partitionEdge*, Less_partitionEdge> pedges; - for(GModel::fiter it = model->firstFace(); it != model->lastFace(); ++it){ - fillit_ ( edgeToElement,(*it)->triangles.begin(),(*it)->triangles.end()); - fillit_ ( edgeToElement,(*it)->quadrangles.begin(),(*it)->quadrangles.end()); - fillit_ ( edgeToElement,(*it)->polygons.begin(),(*it)->polygons.end()); + std::set<partitionVertex*, Less_partitionVertex> pvertices; + std::set<partitionFace*, Less_partitionFace> pfaces; + + // assign partition faces + { + std::multimap<MFace,MElement*,Less_Face> faceToElement; + if (meshDim == 3){ + 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()); + fillit_ ( faceToElement,(*it)->prisms.begin(),(*it)->prisms.end()); + fillit_ ( faceToElement,(*it)->pyramids.begin(),(*it)->pyramids.end()); + fillit_ ( faceToElement,(*it)->polyhedra.begin(),(*it)->polyhedra.end()); + } + } + + // printf("%d edges to elements\n",edgeToElement.size()); + + { + std::multimap<MFace,MElement*,Less_Face>::iterator it = faceToElement.begin(); + Equal_Face oper; + while ( it != faceToElement.end() ){ + MFace e = it->first; + std::vector<MElement*> voe; + do { + voe.push_back(it->second); + ++it; + if (it == faceToElement.end() ) break; + }while (oper (e,it->first)); + assignPartitionBoundary (model,e,pfaces,voe); + } + } + } + + // assign partition edges + { + 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()); + fillit_ ( edgeToElement,(*it)->quadrangles.begin(),(*it)->quadrangles.end()); + fillit_ ( edgeToElement,(*it)->polygons.begin(),(*it)->polygons.end()); + } + } + else if (meshDim == 3){ + for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it){ + fillit_ ( edgeToElement,(*it)->tetrahedra.begin(),(*it)->tetrahedra.end()); + fillit_ ( edgeToElement,(*it)->hexahedra.begin(),(*it)->hexahedra.end()); + fillit_ ( edgeToElement,(*it)->prisms.begin(),(*it)->prisms.end()); + fillit_ ( edgeToElement,(*it)->pyramids.begin(),(*it)->pyramids.end()); + fillit_ ( edgeToElement,(*it)->polyhedra.begin(),(*it)->polyhedra.end()); + } + } + + // printf("%d edges to elements\n",edgeToElement.size()); + + { + std::multimap<MEdge,MElement*,Less_Edge>::iterator it = edgeToElement.begin(); + Equal_Edge oper; + while ( it != edgeToElement.end() ){ + MEdge e = it->first; + std::vector<MElement*> voe; + do { + voe.push_back(it->second); + ++it; + if (it == edgeToElement.end() ) break; + }while (oper (e,it->first)); + assignPartitionBoundary (model,e,pedges,voe,pfaces); + } + } } + // make partition vertices { - std::multimap<MEdge,MElement*,Less_Edge>::iterator it = edgeToElement.begin(); - Equal_Edge oper; - while ( it != edgeToElement.end() ){ - MEdge e = it->first; - std::vector<MElement*> voe; - do { - voe.push_back(it->second); - ++it; - }while (! oper (e,it->first)); - assignPartitionBoundary (model,e,pedges,voe); + 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()); + fillit_ ( vertexToElement,(*it)->quadrangles.begin(),(*it)->quadrangles.end()); + fillit_ ( vertexToElement,(*it)->polygons.begin(),(*it)->polygons.end()); + } + } + else if (meshDim == 3){ + for(GModel::riter it = model->firstRegion(); it != model->lastRegion(); ++it){ + fillit_ ( vertexToElement,(*it)->tetrahedra.begin(),(*it)->tetrahedra.end()); + fillit_ ( vertexToElement,(*it)->hexahedra.begin(),(*it)->hexahedra.end()); + fillit_ ( vertexToElement,(*it)->prisms.begin(),(*it)->prisms.end()); + fillit_ ( vertexToElement,(*it)->pyramids.begin(),(*it)->pyramids.end()); + fillit_ ( vertexToElement,(*it)->polyhedra.begin(),(*it)->polyhedra.end()); + } + } + { + std::multimap<MVertex*,MElement*>::iterator it = vertexToElement.begin(); + while ( it != vertexToElement.end() ){ + MVertex *v = it->first; + std::vector<MElement*> voe; + do { + voe.push_back(it->second); + ++it; + if (it == vertexToElement.end() ) break; + }while (v == it->first); + assignPartitionBoundary (model,v,pvertices,voe,pedges,pfaces); + } } } } - + /******************************************************************************* diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h index a290c602457360f27992a49eda3382674e3883a4..e1e2eb8ce7c880e570604b4b02c817636c4b0213 100644 --- a/Mesh/meshPartitionOptions.h +++ b/Mesh/meshPartitionOptions.h @@ -99,7 +99,7 @@ struct meshPartitionOptions algorithm = 1; edge_matching = 3; refine_algorithm = 3; - createPartitionBoundaries = false;//true; + createPartitionBoundaries = true; } };