From 1772e01edb7d428190a8e13c18caed774b534161 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Mon, 26 Aug 2024 16:37:08 +0200
Subject: [PATCH] update IdentifyPeriodicBoundaries to also support 2d meshes
 (and more stable)

---
 libsrc/meshing/meshclass.cpp   | 46 +++++++++++++++-------------------
 libsrc/meshing/meshclass.hpp   |  4 +--
 libsrc/meshing/python_mesh.cpp |  3 ++-
 3 files changed, 24 insertions(+), 29 deletions(-)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index e324ebba..0727162d 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6812,12 +6812,12 @@ namespace netgen
   //   }
   // #endif
 
-  int Mesh::IdentifyPeriodicBoundaries(const string &s1,
-                                       const string &s2,
+  int Mesh::IdentifyPeriodicBoundaries(const string& id_name,
+                                       const string &s1,
                                        const Transformation<3> &mapping,
                                        double pointTolerance)
   {
-    auto nr = ident->GetMaxNr() + 1;
+    auto nr = ident->GetNr(id_name);
     ident->SetType(nr, Identifications::PERIODIC);
     double lami[4];
     set<int> identified_points;
@@ -6827,44 +6827,38 @@ namespace netgen
         GetBox(pmin, pmax);
         pointTolerance = 1e-8 * (pmax-pmin).Length();
       }
-    for(const auto& se : surfelements)
+    size_t nse = GetDimension() == 3 ? surfelements.Size() : segments.Size();
+    for(auto sei : Range(nse))
       {
-        if(GetBCName(se.index-1) != s1)
+        auto name = GetDimension() == 3 ? GetBCName(surfelements[sei].index-1) :
+          GetBCName(segments[sei].edgenr-1);
+        if(name != s1)
           continue;
 
-        for(const auto& pi : se.PNums())
+        const auto& pnums = GetDimension() == 3 ? surfelements[sei].PNums() :
+          segments[sei].PNums();
+        for(const auto& pi : pnums)
           {
             if(identified_points.find(pi) != identified_points.end())
               continue;
             auto pt = (*this)[pi];
             auto mapped_pt = mapping(pt);
-            // auto other_nr = GetElementOfPoint(mapped_pt, lami, true);
-            auto other_nr = GetSurfaceElementOfPoint(mapped_pt, lami);
-            int index = -1;
-            if(other_nr != 0)
+            bool found = false;
+            for(auto other_pi : Range(points))
               {
-                auto other_el = SurfaceElement(other_nr);
-                for(auto i : Range(other_el.PNums().Size()))
-                  if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < pointTolerance)
-                    {
-                      index = i;
-                      break;
-                    }
-                if(index == -1)
+                if((mapped_pt - (*this)[other_pi]).Length() < pointTolerance)
                   {
-                    cout << "point coordinates = " << pt << endl;
-                    cout << "mapped coordinates = " << mapped_pt << endl;
-                    throw Exception("Did not find mapped point with nr " + ToString(pi) + ", are you sure your mesh is periodic?");
+                    identified_points.insert(pi);
+                    ident->Add(pi, other_pi, nr);
+                    found = true;
+                    break;
                   }
-                auto other_pi = other_el.PNums()[index];
-                identified_points.insert(pi);
-                ident->Add(pi, other_pi, nr);
               }
-            else
+            if(!found)
               {
                 cout << "point coordinates = " << pt << endl;
                 cout << "mapped coordinates = " << mapped_pt << endl;
-                throw Exception("Mapped point with nr " + ToString(pi) + " is outside of mesh, are you sure your mesh is periodic?");
+                throw Exception("Did not find mapped point with nr " + ToString(pi) + ", are you sure your mesh is periodic?");
               }
           }
       }
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 0c7fd57c..41b4832b 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -776,8 +776,8 @@ namespace netgen
     { return facedecoding[i-1]; }      
     // { return facedecoding.Elem(i); }
 
-    int IdentifyPeriodicBoundaries(const string& s1,
-                                   const string& s2,
+    int IdentifyPeriodicBoundaries(const string& id_name,
+                                   const string& s1,
                                    const Transformation<3>& mapping,
                                    double pointTolerance);
 
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 70b762e8..46bd540a 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1247,7 +1247,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("identnr"),
            py::arg("type")=Identifications::PERIODIC)
     .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
-         py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.)
+         py::arg("identification_name"), py::arg("face1"), py::arg("mapping"),
+py::arg("point_tolerance") = -1.)
     .def("GetNrIdentifications", [](Mesh& self)
                                  {
                                    return self.GetIdentifications().GetMaxNr();