[Gmsh] ordering of lines/surfaces after stl merge

[Martin Genet]

On Tue 03 Sep 2013 01:42:59 AM PDT, Christophe Geuzaine wrote:
> On 02 Sep 2013, at 19:38, Martin Genet <martin.genet at ucsfmedicalcenter.org> wrote:
>> On 08/31/2013 01:19 AM, Christophe Geuzaine wrote:
>>> On 29 Aug 2013, at 20:44, Martin Genet <martin.genet at ucsfmedicalcenter.org> wrote:
>>>> Sorry to insist, but I'm really stuck here. Hopefully someone can help me understand what's going on.
>>>>
>>>> I made a minimal example (minimal.geo), to illustrate my problem. I first merge the a couple STL files (STL.tar.gz), then create a couple Compound Lines, and then eventually a Line Loop. However, the order of Compound Lines seems to depend on whether or not I create the Line Loop (see lines-WithoutLineLoop & lines-WithLineLoop). Is that the expected behavior? Is there a logic behind this?
>>>
>>> Hi Martin - It's a limitation of the current implementation of "CreateTopology". In this routine (which reconstructs all the connectivity data that is missing in the STL), some operations depend on sets of entities (vertices, elements), which are sorted by their memory address (pointer). The order in which features (curves) are created depends on this ordering, which itself changes from run to run.
>>>
>>> We should definitely change this, as it leads to numbering are change from run to run. It's in the TODO list...
>>>
>>> Christophe
>>
>> Thanks Christophe. Is it something a non-GMSH expert could work out?
>> I'll give it a try. Do you have any instructions? That would help a lot…
>
> The code is quite ugly, but might still be readable :-) It's in gmsh/Geo/GModel.cpp. The routine to look at is GModel::createTopologyFromMesh(), as well as all the subroutines called from there.

I've tried to add vectors along with sets for faces and edges in GModel, in order to keep track of the order in which they are added, and then iterate over these vectors instead of the sets to create the entities. I'm attaching a small patch, showing this. However, it doesn't seem to have any effect on my problem, so I guess I didn't correct the right functions. Let me know if you have additional instructions regarding this issue, otherwise I'll wait for the fix. Thanks again so much for your help.

Martin

-------------- next part --------------
Index: Geo/GModel.h
===================================================================
--- Geo/GModel.h	(revision 16551)
+++ Geo/GModel.h	(working copy)
@@ -150,6 +150,10 @@
   std::set<GEdge*, GEntityLessThan> edges;
   std::set<GVertex*, GEntityLessThan> vertices;
 
+  // the vectors of faces and edges in the model (ordered by insertion time)
+  std::vector<GFace*> ordered_faces;
+  std::vector<GEdge*> ordered_edges;
+
   // map between the pair <dimension, elementary or physical number>
   // and an optional associated name
   std::map<std::pair<int, int>, std::string> physicalNames, elementaryNames;
@@ -266,6 +270,16 @@
   eiter lastEdge() { return edges.end(); }
   viter lastVertex() { return vertices.end(); }
 
+  // ordered face and edge iterators
+  typedef std::vector<GFace*>::iterator ofiter;
+  typedef std::vector<GEdge*>::iterator oeiter;
+
+  // get an iterator initialized to the first/last ordered entity in this model
+  ofiter firstOrderedFace() { return ordered_faces.begin(); }
+  oeiter firstOrderedEdge() { return ordered_edges.begin(); }
+  ofiter lastOrderedFace() { return ordered_faces.end(); }
+  oeiter lastOrderedEdge() { return ordered_edges.end(); }
+
   // find the entity with the given tag
   GRegion *getRegionByTag(int n) const;
   GFace *getFaceByTag(int n) const;
@@ -275,8 +289,8 @@
 
   // add/remove an entity in the model
   void add(GRegion *r) { regions.insert(r); }
-  void add(GFace *f) { faces.insert(f); }
-  void add(GEdge *e) { edges.insert(e); }
+  void add(GFace *f) { std::pair<fiter, bool> res; res = faces.insert(f); ordered_faces.push_back(*res.first);}
+  void add(GEdge *e) { std::pair<eiter, bool> res; res = edges.insert(e); ordered_edges.push_back(*res.first);}
   void add(GVertex *v) { vertices.insert(v); }
   void remove(GRegion *r);
   void remove(GFace *f);
Index: Geo/GModelIO_GEO.cpp
===================================================================
--- Geo/GModelIO_GEO.cpp	(revision 16551)
+++ Geo/GModelIO_GEO.cpp	(working copy)
@@ -57,7 +57,7 @@
     Tree_Add(this->getGEOInternals()->Points, &v);
   }
 
-  for(eiter it = firstEdge(); it != lastEdge(); it++){
+  for(oeiter it = firstOrderedEdge(); it != lastOrderedEdge(); it++){
     if((*it)->geomType() == GEntity::DiscreteCurve){
       Curve *c = Create_Curve((*it)->tag(), MSH_SEGM_DISCRETE, 1,
                               NULL, NULL, -1, -1, 0., 1.);
@@ -85,7 +85,7 @@
     }
   }
 
-  for(fiter it = firstFace(); it != lastFace(); it++){
+  for(ofiter it = firstOrderedFace(); it != lastOrderedFace(); it++){
     if((*it)->geomType() == GEntity::DiscreteSurface){
       Surface *s = Create_Surface((*it)->tag(), MSH_SURF_DISCRETE);
       std::list<GEdge*> edges = (*it)->edges();