new algo to compute partition boundaries (much faster)

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;