Merge branch 'get_sub_mesh' into 'master'

Utility function to extract a subregion of the mesh

See merge request ngsolve/netgen!644

libsrc/meshing/meshclass.cpp

@@ -1,5 +1,6 @@
 #include <mystdlib.h>
 #include <atomic>
+#include <regex>
 #include <set>
 #include "meshing.hpp"
 #include "../general/gzstream.h"
@@ -319,26 +320,27 @@ namespace netgen
     hmin = mesh2.hmin;
     maxhdomain = mesh2.maxhdomain;
 
+    // Remap string* values to new mesh
+    std::map<const string*, string*> names_map;
+    for (auto fi : Range(facedecoding))
+      names_map[&mesh2.facedecoding[fi].bcname] = &facedecoding[fi].bcname;
 
     materials.SetSize( mesh2.materials.Size() );
     for ( int i = 0; i < mesh2.materials.Size(); i++ )
-      if ( mesh2.materials[i] ) materials[i] = new string ( *mesh2.materials[i] );
+    {
+      const string * old_name = mesh2.materials[i];
+      if ( old_name ) materials[i] = dimension == 2 ? names_map[old_name] : new string ( *old_name );
       else materials[i] = 0;
+    }
 
-    std::map<const string*, string*> bcmap;
     bcnames.SetSize( mesh2.bcnames.Size() );
     for ( int i = 0; i < mesh2.bcnames.Size(); i++ )
     {
-      if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] );
+      const string * old_name = mesh2.bcnames[i];
+      if ( old_name ) bcnames[i] = dimension == 3 ? names_map[old_name] : new string ( *old_name );
       else bcnames[i] = 0;
-      bcmap[mesh2.bcnames[i]] = bcnames[i];
     }
 
-    // Remap string* members in FaceDescriptor to new mesh
-    for (auto & f : facedecoding)
-      f.SetBCName( bcmap[&f.GetBCName()] );
-
-
     cd2names.SetSize(mesh2.cd2names.Size());
     for (int i=0; i < mesh2.cd2names.Size(); i++)
       if (mesh2.cd2names[i]) cd2names[i] = new string(*mesh2.cd2names[i]);
@@ -7199,6 +7201,83 @@ namespace netgen
     UpdateTopology();
   }
 
+  shared_ptr<Mesh> Mesh :: GetSubMesh(string domains, string faces) const
+  {
+    // Copy the mesh into a new one, then delete unwanted elements
+    // Unused points are deleted by the Compress() function at the end
+    auto mesh_ptr = make_unique<Mesh>();
+    auto & mesh = *mesh_ptr;
+    mesh = (*this);
+
+    auto ndomains = GetNDomains();
+    auto nfaces = GetNFD();
+
+    BitArray keep_point(GetNP()+1);
+    BitArray keep_face(nfaces+1);
+    BitArray keep_domain(ndomains+1);
+    keep_point.Clear();
+    keep_face.Clear();
+    keep_domain.Clear();
+
+    regex regex_faces(faces);
+    regex regex_domains(domains);
+
+    for(auto dom : Range(ndomains))
+    {
+      if(regex_match(mesh.GetMaterial(dom), regex_domains))
+        keep_domain.SetBit(dom);
+    }
+
+    for(auto fi : Range(nfaces))
+    {
+      auto & fd = mesh.FaceDescriptors()[fi];
+      if (keep_domain[fd.DomainIn()] || keep_domain[fd.DomainOut()] || regex_match(fd.GetBCName(), regex_faces))
+        keep_face.SetBit(fd.BCProperty());
+    }
+
+    auto filter_elements = [&mesh, &keep_point](auto & elements, auto & keep_region)
+    {
+      for(auto & el : elements)
+      {
+        if(keep_region[el.GetIndex()])
+          for (auto pi : el.PNums())
+            keep_point.SetBit(pi);
+        else
+          el.Delete();
+      }
+    };
+
+    filter_elements(mesh.VolumeElements(), keep_domain);
+    filter_elements(mesh.SurfaceElements(), keep_face);
+
+    // Keep line segments only if all points are kept
+    // Check them in reverse order because they are deleted from the end
+    auto nsegments = mesh.LineSegments().Size();
+    for(auto i : Range(nsegments))
+    {
+      SegmentIndex segi = nsegments-i-1;
+      auto seg = mesh[segi];
+      bool keep = true;
+      for(auto pi : seg.PNums())
+        keep &= keep_point[pi];
+
+      if(!keep)
+        mesh.LineSegments().DeleteElement(segi);
+    }
+
+    // Check in reverse order because they are deleted from the end
+    auto npointelements = mesh.pointelements.Size();
+    for(auto i : Range(npointelements))
+    {
+      auto pel = mesh.pointelements[npointelements-i-1];
+      if(!keep_point[pel.pnum])
+        mesh.pointelements.DeleteElement(npointelements-i-1);
+    }
+
+    mesh.Compress();
+    return mesh_ptr;
+  }
+
   void Mesh :: SetMaterial (int domnr, const string & mat)
   {
     if (domnr > materials.Size())

libsrc/meshing/meshclass.hpp

@@ -700,6 +700,7 @@ namespace netgen
     void SetCommunicator(NgMPI_Comm acomm);
     
     DLL_HEADER void SplitFacesByAdjacentDomains();
+    DLL_HEADER shared_ptr<Mesh> GetSubMesh(string domains="", string faces="") const;
 
     ///
     DLL_HEADER void SetMaterial (int domnr, const string & mat);

libsrc/meshing/python_mesh.cpp

@@ -1251,6 +1251,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def ("GetCD3Name", &Mesh::GetCD3Name)
     .def ("SetCD3Name", &Mesh::SetCD3Name)
     .def ("SplitFacesByAdjacentDomains", &Mesh::SplitFacesByAdjacentDomains)
+    .def ("GetSubMesh", &Mesh::GetSubMesh, py::arg("domains")="", py::arg("faces")="")
     .def("GetIdentifications", [](Mesh & self) -> py::list
          {
            py::list points;