diff --git a/Plugin/SimplePartition.cpp b/Plugin/SimplePartition.cpp
index 9495e425e748ee4cc11f5e4a0196b4d3128ee55a..62988f6463c20d03a081abf66a72a9714ba68f36 100644
--- a/Plugin/SimplePartition.cpp
+++ b/Plugin/SimplePartition.cpp
@@ -5,7 +5,10 @@
#include "SimplePartition.h"
#include "GModel.h"
+#include "partitionFace.h"
+#include "partitionEdge.h"
#include "MElement.h"
+#include "MLine.h"
#include "MFace.h"
#include "MEdge.h"
#if defined(HAVE_MESH)
@@ -57,34 +60,40 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
GModel *m = GModel::current();
SBoundingBox3d bbox = m->bounds();
double pmin = bbox.min()[direction], pmax = bbox.max()[direction];
+ std::vector<double> pp(numSlices + 1);
+ for(int p = 0; p <= numSlices; p++)
+ pp[p] = pmin + p * (pmax - pmin) / numSlices;
int dim = m->getDim();
std::vector<GEntity*> entities;
m->getEntities(entities);
- std::map<int, std::set<MFace, Less_Face> > bndFaces;
- std::map<int, std::set<MEdge, Less_Edge> > bndEdges;
+ std::vector<std::set<MFace, Less_Face> > bndFaces(numSlices);
+ std::vector<std::set<MEdge, Less_Edge> > bndEdges(numSlices);
+ std::vector<std::set<MElement*> > cutElements(numSlices);
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
if(ge->dim() != dim) continue;
for(int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
- double val = pmax;
+ double valmin = pmax;
+ double valmax = pmin;
bool bnd = false;
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
- val = std::min(val, v->point()[direction]);
+ valmin = std::min(valmin, v->point()[direction]);
+ valmax = std::max(valmax, v->point()[direction]);
if(v->onWhat() && v->onWhat()->dim() != dim) bnd = true;
}
for(int p = 0; p < numSlices; p++){
- double p1 = pmin + p * (pmax - pmin) / numSlices;
- double p2 = pmin + (p + 1) * (pmax - pmin) / numSlices;
- if(val >= p1 && val < p2){
+ if(valmin >= pp[p] && valmin < pp[p + 1]){
e->setPartition(p + 1);
if(bnd){
for(int k = 0; k < e->getNumEdges(); k++)
- bndEdges[p + 1].insert(e->getEdge(k));
+ bndEdges[p].insert(e->getEdge(k));
for(int k = 0; k < e->getNumFaces(); k++)
- bndFaces[p + 1].insert(e->getFace(k));
+ bndFaces[p].insert(e->getFace(k));
}
+ if(valmax >= pp[p + 1])
+ cutElements[p].insert(e);
break;
}
}
@@ -94,6 +103,8 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
// partition lower dimension elements
Msg::Info("Partitioning lower dimension elements");
+ std::set<MFace, Less_Face> allLowerDimFaces;
+ std::set<MEdge, Less_Edge> allLowerDimEdges;
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *ge = entities[i];
if(ge->dim() == dim) continue;
@@ -101,8 +112,9 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
MElement *e = ge->getMeshElement(j);
if(e->getDim() == 2){
MFace f = e->getFace(0);
+ if(createBoundaries) allLowerDimFaces.insert(f);
for(int p = 0; p < numSlices; p++){
- if(bndFaces[p + 1].find(f) != bndFaces[p + 1].end()){
+ if(bndFaces[p].find(f) != bndFaces[p].end()){
e->setPartition(p + 1);
break;
}
@@ -110,8 +122,9 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
}
else if(e->getDim() == 1){
MEdge f = e->getEdge(0);
+ if(createBoundaries) allLowerDimEdges.insert(f);
for(int p = 0; p < numSlices; p++){
- if(bndEdges[p + 1].find(f) != bndEdges[p + 1].end()){
+ if(bndEdges[p].find(f) != bndEdges[p].end()){
e->setPartition(p + 1);
break;
}
@@ -121,13 +134,9 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
}
if(createBoundaries){
-#if defined(HAVE_MESH)
Msg::Info("Creating partition boundaries");
+#if 0 && defined(HAVE_MESH) // correct & general, but too slow
CreatePartitionBoundaries(m, false, false);
-#else
- Msg::Error("Creating partition boundaries requires the mesh module");
-#endif
- // renumber partition boundaries using the natural slicing order
Msg::Info("Renumbering partition boundaries");
std::vector<std::pair<double, GFace*> > parFaces;
std::vector<std::pair<double, GEdge*> > parEdges;
@@ -136,22 +145,22 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
GEntity *ge = entities[i];
if(ge->geomType() == GEntity::PartitionSurface ||
ge->geomType() == GEntity::PartitionCurve){
- double val = pmax;
+ double valmin = pmax;
for(int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
- val = std::min(val, v->point()[direction]);
+ valmin = std::min(valmin, v->point()[direction]);
}
}
if(ge->dim() == 2){
GFace *gf = (GFace*)ge;
- parFaces.push_back(std::pair<double, GFace*>(val, gf));
+ parFaces.push_back(std::pair<double, GFace*>(valmin, gf));
m->remove(gf);
}
else{
GEdge *gc = (GEdge*)ge;
- parEdges.push_back(std::pair<double, GEdge*>(val, gc));
+ parEdges.push_back(std::pair<double, GEdge*>(valmin, gc));
m->remove(gc);
}
}
@@ -170,6 +179,87 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
ge->setMeshMaster(-i-1);
m->add(ge);
}
+#else
+ for(int p = 0; p < numSlices - 1; p++){
+ std::vector<int> v2(2);
+ v2[0] = p + 1;
+ v2[1] = p + 2;
+ if(dim == 2){
+ // create partition edge
+ partitionEdge *pe = new partitionEdge(m, -(p + 1), 0, 0, v2);
+ m->add(pe);
+ // compute boundary of cut surface elements
+ std::set<MEdge, Less_Edge> myBndEdges;
+ for(std::set<MElement*>::iterator it = cutElements[p].begin();
+ it != cutElements[p].end(); it++){
+ for(int i = 0; i < (*it)->getNumEdges(); i++){
+ MEdge e = (*it)->getEdge(i);
+ if(myBndEdges.find(e) == myBndEdges.end())
+ myBndEdges.insert(e);
+ else
+ myBndEdges.erase(e);
+ }
+ }
+ // keep edges whose vertices are all >= than the plane, but are not on a
+ // curve of the original mesh
+ for(std::set<MEdge, Less_Edge>::iterator it = myBndEdges.begin();
+ it != myBndEdges.end(); it++){
+ bool bnd = true;
+ for(int j = 0; j < it->getNumVertices(); j++){
+ if(it->getVertex(j)->point()[direction] < pp[p + 1]){
+ bnd = false;
+ break;
+ }
+ }
+ if(bnd){
+ if(allLowerDimEdges.find(*it) == allLowerDimEdges.end()){
+ pe->lines.push_back(new MLine(it->getVertex(0), it->getVertex(1)));
+ }
+ }
+ }
+ }
+ if(dim == 3){
+ // create partition face
+ partitionFace *pf = new partitionFace(m, -(p + 1), v2);
+ m->add(pf);
+ // compute boundary of cut elements
+ std::set<MFace, Less_Face> myBndFaces;
+ for(std::set<MElement*>::iterator it = cutElements[p].begin();
+ it != cutElements[p].end(); it++){
+ for(int i = 0; i < (*it)->getNumFaces(); i++){
+ MFace f = (*it)->getFace(i);
+ if(myBndFaces.find(f) == myBndFaces.end())
+ myBndFaces.insert(f);
+ else
+ myBndFaces.erase(f);
+ }
+ }
+ // keep faces whose vertices are all >= than the plane, but are not on a
+ // surface of the original mesh
+ for(std::set<MFace, Less_Face>::iterator it = myBndFaces.begin();
+ it != myBndFaces.end(); it++){
+ bool bnd = true;
+ for(int j = 0; j < it->getNumVertices(); j++){
+ if(it->getVertex(j)->point()[direction] < pp[p + 1]){
+ bnd = false;
+ break;
+ }
+ }
+ if(bnd){
+ if(allLowerDimFaces.find(*it) == allLowerDimFaces.end()){
+ if (it->getNumVertices() == 3)
+ pf->triangles.push_back
+ (new MTriangle(it->getVertex(0), it->getVertex(1), it->getVertex(2)));
+ else
+ pf->quadrangles.push_back
+ (new MQuadrangle(it->getVertex(0), it->getVertex(1),
+ it->getVertex(2), it->getVertex(3)));
+ }
+ }
+ }
+ }
+ }
+#endif
}
return v;