better: correct lower dim + sort partition boundaries w.r.t. slices

Geo/GEntity.cpp

@@ -72,6 +72,8 @@ GRegion *GEntity::cast2Region() { return dynamic_cast<GRegion*>(this); }
 // sets the entity m from which the mesh will be copied
 void GEntity::setMeshMaster(int m_signed){
 
+  if(m_signed == tag()){ _meshMaster = m_signed; return; }
+
   //  printf("setting mesh master %d to mesh entity %d\n",m_signed,tag());
 
   GEntity *gMaster = 0;
@@ -81,12 +83,12 @@ void GEntity::setMeshMaster(int m_signed){
   case 1 : gMaster = model()->getEdgeByTag(m); break;
   case 2 : gMaster = model()->getFaceByTag(m); break;
   case 3 : gMaster = model()->getRegionByTag(m); break;
-  } 
+  }
   if (!gMaster){
     Msg::Fatal("Model entity %d of dimension %d cannot be the mesh master of model entity %d",m,dim(), tag());
   }
   int masterOfMaster = gMaster->meshMaster();
-  
+
   if (masterOfMaster == gMaster->tag()){
     _meshMaster = m_signed;
   }
@@ -98,7 +100,7 @@ void GEntity::setMeshMaster(int m_signed){
 
 // gets the entity from which the mesh will be copied
 int GEntity::meshMaster() const{
-  
+
   if (_meshMaster == tag()) return tag();
 
   GEntity *gMaster = 0;
@@ -107,12 +109,12 @@ int GEntity::meshMaster() const{
   case 1 : gMaster = model()->getEdgeByTag(abs(_meshMaster)); break;
   case 2 : gMaster = model()->getFaceByTag(abs(_meshMaster)); break;
   case 3 : gMaster = model()->getRegionByTag(abs(_meshMaster)); break;
-  } 
+  }
   if (!gMaster){
     Msg::Fatal("meshMaster : Model entity %d of dimension %d cannot be the mesh master of model entity %d",_meshMaster,dim(),tag());
   }
   int masterOfMaster = gMaster->meshMaster();
-  
+
   if (masterOfMaster == gMaster->tag()){
     return _meshMaster ;
   }

Plugin/SimplePartition.cpp

@@ -6,6 +6,8 @@
 #include "SimplePartition.h"
 #include "GModel.h"
 #include "MElement.h"
+#include "MFace.h"
+#include "MEdge.h"
 #if defined(HAVE_MESH)
 #include "meshPartition.h"
 #endif
@@ -45,26 +47,40 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
   int numSlices = (int)SimplePartitionOptions_Number[0].def;
   int direction = (int)SimplePartitionOptions_Number[1].def;
 
-  // partition the current model
+  // partition the highest dimension elements in the current model (lower
+  // dimension elements en boundaries cannot be tagged a priori: there are
+  // special geometrical cases)
   GModel *m = GModel::current();
   SBoundingBox3d bbox = m->bounds();
   double pmin = bbox.min()[direction], pmax = bbox.max()[direction];
+  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;
   for(unsigned int i = 0; i < entities.size(); i++){
-    int n = entities[i]->getNumMeshElements();
-    for(int j = 0; j < n; j++){
-      MElement *e = entities[i]->getMeshElement(j);
+    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;
+      bool bnd = false;
       for(int k = 0; k < e->getNumVertices(); k++){
         MVertex *v = e->getVertex(k);
         val = std::min(val, 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){
           e->setPartition(p + 1);
+          if(bnd){
+            for(int k = 0; k < e->getNumEdges(); k++)
+              bndEdges[p + 1].insert(e->getEdge(k));
+            for(int k = 0; k < e->getNumFaces(); k++)
+              bndFaces[p + 1].insert(e->getFace(k));
+          }
           break;
         }
       }
@@ -72,8 +88,82 @@ PView *GMSH_SimplePartitionPlugin::execute(PView *v)
   }
   m->recomputeMeshPartitions();
 
+  // partition lower dimension elements
+  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);
+      if(e->getDim() == 2){
+        MFace f = e->getFace(0);
+        for(int p = 0; p < numSlices; p++){
+          if(bndFaces[p + 1].find(f) != bndFaces[p + 1].end()){
+            e->setPartition(p + 1);
+            break;
+          }
+        }
+      }
+      else if(e->getDim() == 1){
+        MEdge f = e->getEdge(0);
+        for(int p = 0; p < numSlices; p++){
+          if(bndEdges[p + 1].find(f) != bndEdges[p + 1].end()){
+            e->setPartition(p + 1);
+            break;
+          }
+        }
+      }
+    }
+  }
+
+  // create partition boundaries
 #if defined(HAVE_MESH)
   CreatePartitionBoundaries(m, false, false);
 #endif
+
+  // renumber partition boundaries using the natural slicing order
+  std::vector<std::pair<double, GFace*> > parFaces;
+  std::vector<std::pair<double, GEdge*> > parEdges;
+  m->getEntities(entities);
+  for(unsigned int i = 0; i < entities.size(); i++){
+    GEntity *ge = entities[i];
+    if(ge->geomType() == GEntity::PartitionSurface ||
+       ge->geomType() == GEntity::PartitionCurve){
+      double val = 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]);
+        }
+      }
+      if(ge->dim() == 2){
+        GFace *gf = (GFace*)ge;
+        parFaces.push_back(std::pair<double, GFace*>(val, gf));
+        m->remove(gf);
+      }
+      else{
+        GEdge *gc = (GEdge*)ge;
+        parEdges.push_back(std::pair<double, GEdge*>(val, gc));
+        m->remove(gc);
+      }
+    }
+  }
+
+  std::sort(parFaces.begin(), parFaces.end());
+  for(unsigned int i = 0; i < parFaces.size(); i++){
+    GFace *gf = parFaces[i].second;
+    gf->setTag(-i-1);
+    gf->setMeshMaster(-i-1);
+    m->add(gf);
+  }
+
+  std::sort(parEdges.begin(), parEdges.end());
+  for(unsigned int i = 0; i < parEdges.size(); i++){
+    GEdge *ge = parEdges[i].second;
+    ge->setTag(-i-1);
+    ge->setMeshMaster(-i-1);
+    m->add(ge);
+  }
+
   return v;
 }