From 4e31878f89443dce169638f81af111da44123d35 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Wed, 31 Jan 2024 11:55:53 +0100
Subject: [PATCH] Utility function to split faces when they have more than two
 adjacent domains

---
 libsrc/meshing/meshclass.cpp   | 47 ++++++++++++++++++++++++++++++++++
 libsrc/meshing/meshclass.hpp   |  2 ++
 libsrc/meshing/python_mesh.cpp |  1 +
 3 files changed, 50 insertions(+)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 382b669b..c7b70efc 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -7152,6 +7152,53 @@ namespace netgen
     SetNextMajorTimeStamp();
   }
 
+  void Mesh :: SplitFacesByAdjacentDomains ()
+  {
+    UpdateTopology();
+    std::map<std::tuple<int, int, int>, int> face_doms_2_new_face;
+    int nfaces = FaceDescriptors().Size();
+    Array<bool> first_visit(nfaces);
+    first_visit = true;
+
+    for (auto sei : Range(SurfaceElements()))
+      {
+        int eli0, eli1;
+        GetTopology().GetSurface2VolumeElement(sei+1, eli0, eli1);
+        auto & sel = (*this)[sei];
+        int face = sel.GetIndex();
+        int domin = VolumeElement(eli0).GetIndex();
+        int domout = eli1 ? VolumeElement(eli1).GetIndex() : 0;
+        if(domin < domout)
+          swap(domin, domout);
+        auto key = std::make_tuple(face, domin, domout);
+        if(face_doms_2_new_face.find(key) == face_doms_2_new_face.end())
+          {
+            if(!first_visit[face-1]) {
+              nfaces++;
+              FaceDescriptor new_fd = FaceDescriptors()[face-1];
+              new_fd.bcprop = nfaces;
+              new_fd.domin = domin;
+              new_fd.domout = domout;
+              AddFaceDescriptor(new_fd);
+              SetBCName(nfaces-1, new_fd.GetBCName());
+              face_doms_2_new_face[key] = nfaces;
+            }
+            else {
+              face_doms_2_new_face[key] = face;
+              auto & fd = FaceDescriptors()[face-1];
+              fd.domin = domin;
+              fd.domout = domout;
+            }
+            first_visit[face-1] = false;
+          }
+          sel.SetIndex(face_doms_2_new_face[key]);
+      }
+    SetNextMajorTimeStamp();
+    RebuildSurfaceElementLists ();
+    CalcSurfacesOfNode();
+    UpdateTopology();
+  }
+
   void Mesh :: SetMaterial (int domnr, const string & mat)
   {
     if (domnr > materials.Size())
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index a1adf171..e946244f 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -699,6 +699,8 @@ namespace netgen
     auto & GetCommunicator() const { return this->comm; }
     void SetCommunicator(NgMPI_Comm acomm);
     
+    DLL_HEADER void SplitFacesByAdjacentDomains();
+
     ///
     DLL_HEADER void SetMaterial (int domnr, const string & mat);
     ///
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 910a5e13..71bb1245 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1250,6 +1250,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 
     .def ("GetCD3Name", &Mesh::GetCD3Name)
     .def ("SetCD3Name", &Mesh::SetCD3Name)
+    .def ("SplitFacesByAdjacentDomains", &Mesh::SplitFacesByAdjacentDomains)
     .def("GetIdentifications", [](Mesh & self) -> py::list
          {
            py::list points;