From b8d313f056459e3a0ee4c0c498f50b10b63339e3 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Sun, 15 Mar 2020 18:02:50 +0100
Subject: [PATCH 1/2] identify periodic boundaries

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

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 8c340041..5c774dfa 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -1,5 +1,6 @@
 #include <mystdlib.h>
 #include <atomic>
+#include <set>
 #include "meshing.hpp"
 
 #ifdef NG_PYTHON
@@ -6112,7 +6113,48 @@ namespace netgen
   //   }
   // #endif
 
+  int Mesh::IdentifyPeriodicBoundaries(const string &s1,
+                                       const string &s2,
+                                       const Transformation<3> &mapping)
+  {
+    auto nr = ident->GetMaxNr() + 1;
+    double lami[4];
+    GetElementOfPoint({0,0,0}, lami, true);
+    set<int> identified_points;
+    Point3d pmin, pmax;
+    GetBox(pmin, pmax);
+    auto eps = 1e-10 * (pmax-pmin).Length();
+    for(const auto& se : surfelements)
+      {
+        if(GetBCName(se.index-1) != s1)
+          continue;
 
+        for(const auto& pi : se.PNums())
+          {
+            // cout << "pi = " << pi << endl;
+            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);
+            int index = -1;
+            auto other_el = VolumeElement(other_nr);
+            for(auto i : Range(4))
+              if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < eps)
+                {
+                  index = i;
+                  break;
+                }
+            if(index == -1)
+              throw Exception("Did not find mapped point, are you sure your mesh is periodic?");
+            auto other_pi = other_el.PNums()[index];
+            identified_points.insert(pi);
+            ident->Add(pi, other_pi, nr);
+            // cout << "other pi = " << other_pi << endl;
+          }
+      }
+    return nr;
+  }
 
   void Mesh :: InitPointCurve(double red, double green, double blue) const
   {
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index e455f02e..f312a5da 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -710,6 +710,10 @@ namespace netgen
     FaceDescriptor & GetFaceDescriptor (int i) 
     { return facedecoding.Elem(i); }
 
+    int IdentifyPeriodicBoundaries(const string& s1,
+                                   const string& s2,
+                                   const Transformation<3>& mapping);
+
     // #ifdef NONE
     //   /*
     //     Identify points pi1 and pi2, due to
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 23cb5d3e..b0ac0c36 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -855,6 +855,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("pid2"),
            py::arg("identnr"),
            py::arg("type"))
+    .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries)
     .def ("CalcLocalH", &Mesh::CalcLocalH)
     .def ("SetMaxHDomain", [] (Mesh& self, py::list maxhlist)
           {

From ff60ca3f554cd22e85e274ddf02def6c8249f751 Mon Sep 17 00:00:00 2001
From: Christopher Lackner <clackner@cerbsim.com>
Date: Tue, 17 Mar 2020 15:32:42 +0100
Subject: [PATCH 2/2] fix identify periodic

---
 libsrc/meshing/meshclass.cpp   | 8 +++-----
 libsrc/meshing/python_mesh.cpp | 4 ++++
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 5c774dfa..5441b732 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6118,12 +6118,12 @@ namespace netgen
                                        const Transformation<3> &mapping)
   {
     auto nr = ident->GetMaxNr() + 1;
+    ident->SetType(nr, Identifications::PERIODIC);
     double lami[4];
-    GetElementOfPoint({0,0,0}, lami, true);
     set<int> identified_points;
     Point3d pmin, pmax;
     GetBox(pmin, pmax);
-    auto eps = 1e-10 * (pmax-pmin).Length();
+    auto eps = 1e-8 * (pmax-pmin).Length();
     for(const auto& se : surfelements)
       {
         if(GetBCName(se.index-1) != s1)
@@ -6131,12 +6131,11 @@ namespace netgen
 
         for(const auto& pi : se.PNums())
           {
-            // cout << "pi = " << pi << endl;
             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);
+            auto other_nr = GetElementOfPoint(mapped_pt, lami, true);
             int index = -1;
             auto other_el = VolumeElement(other_nr);
             for(auto i : Range(4))
@@ -6150,7 +6149,6 @@ namespace netgen
             auto other_pi = other_el.PNums()[index];
             identified_points.insert(pi);
             ident->Add(pi, other_pi, nr);
-            // cout << "other pi = " << other_pi << endl;
           }
       }
     return nr;
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index b0ac0c36..e5a8c5cc 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -856,6 +856,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("identnr"),
            py::arg("type"))
     .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries)
+    .def("GetNrIdentifications", [](Mesh& self)
+                                 {
+                                   return self.GetIdentifications().GetMaxNr();
+                                 })
     .def ("CalcLocalH", &Mesh::CalcLocalH)
     .def ("SetMaxHDomain", [] (Mesh& self, py::list maxhlist)
           {