From 0d05063657df452099c84e61aada66e732cb98b0 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 24 Jan 2020 10:26:09 +0100
Subject: [PATCH] move physical group definitions to new discrete entity
 created by compound meshing constraint, except for surfaces when we do the
 "magic" reclassification on the original surfaces

---
 CHANGELOG.txt           | 3 ++-
 Common/DefaultOptions.h | 2 +-
 Geo/GEdge.cpp           | 4 ++++
 Geo/GFace.cpp           | 6 ++++++
 4 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 4b358d26cb..0e8f069825 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,5 +1,6 @@
 4.5.2 (Work-in-progress): periodic meshes now obey reorientation constraints;
-small bug fixes.
+physical group definitions now follow compound meshing constraints; small bug
+fixes.
 
 4.5.1 (December 28, 2019): new Min and Max commands in .geo files;
 Mesh.MinimumCirclePoints now behaves the same with all geometry kernels; fixed
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index ae55785bd7..a531fa1800 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1061,7 +1061,7 @@ StringXNumber MeshOptions_Number[] = {
     "entity, 3: by partition)" },
   { F|O, "CompoundClassify" , opt_mesh_compound_classify , 1. ,
     "How are surface mesh elements classified on compounds? (0: on the new discrete "
-    "entity, 1: on the original geometrical entity - incompatible with e.g. high-order "
+    "surface, 1: on the original geometrical surfaces - incompatible with e.g. high-order "
     "meshing)" },
   { F|O, "CompoundCharacteristicLengthFactor" , opt_mesh_compound_lc_factor , 0.5 ,
     "Mesh size factor applied to compound parts" },
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 3e102edaea..9a2dc11a8b 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -738,6 +738,7 @@ static void meshCompound(GEdge *ge)
   // no new mesh nodes are created here
   discreteEdge *de = new discreteEdge(ge->model(), ge->tag() + 100000);
   ge->model()->add(de);
+  std::vector<int> phys;
   for(std::size_t i = 0; i < ge->compound.size(); i++) {
     GEdge *c = (GEdge *)ge->compound[i];
     // cannot use the same line elements, as they get deleted in createGeometry
@@ -746,6 +747,8 @@ static void meshCompound(GEdge *ge)
                                     c->lines[j]->getVertex(1)));
     }
     c->compoundCurve = de;
+    phys.insert(phys.end(), c->physicals.begin(), c->physicals.end());
+    c->physicals.clear();
   }
   // create the geometry of the compound
   de->createGeometry(true);
@@ -758,6 +761,7 @@ static void meshCompound(GEdge *ge)
   de->deleteVertexArrays();
   // mesh the compound
   de->mesh(false);
+  de->physicals = phys;
 }
 #endif
 
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 2f4622d920..cac83eb3d7 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1486,6 +1486,7 @@ static void meshCompound(GFace *gf, bool verbose)
 
   std::set<GEdge*, GEntityPtrLessThan> bnd, emb1;
   std::set<GVertex*, GEntityPtrLessThan> emb0;
+  std::vector<int> phys;
   for(std::size_t i = 0; i < gf->compound.size(); i++) {
     GFace *c = (GFace *)gf->compound[i];
     df->triangles.insert(df->triangles.end(), c->triangles.begin(),
@@ -1515,6 +1516,10 @@ static void meshCompound(GFace *gf, bool verbose)
       c->mesh_vertices.clear();
     }
     c->compoundSurface = df;
+    if(!magic) {
+      phys.insert(phys.end(), c->physicals.begin(), c->physicals.end());
+      c->physicals.clear();
+    }
   }
 
   std::set<GEdge*, GEntityPtrLessThan> bndc;
@@ -1550,6 +1555,7 @@ static void meshCompound(GFace *gf, bool verbose)
   df->mesh(verbose);
 
   if(!magic){
+    df->physicals = phys;
     return;
   }
 
-- 
GitLab