Merge branch 'identify_periodic' into 'master'

Identify periodic

See merge request jschoeberl/netgen!316

libsrc/meshing/meshclass.cpp

@@ -1,5 +1,6 @@
 #include <mystdlib.h>
 #include <atomic>
+#include <set>
 #include "meshing.hpp"
 
 #ifdef NG_PYTHON
@@ -6112,7 +6113,46 @@ namespace netgen
   //   }
   // #endif
 
+  int Mesh::IdentifyPeriodicBoundaries(const string &s1,
+                                       const string &s2,
+                                       const Transformation<3> &mapping)
+  {
+    auto nr = ident->GetMaxNr() + 1;
+    ident->SetType(nr, Identifications::PERIODIC);
+    double lami[4];
+    set<int> identified_points;
+    Point3d pmin, pmax;
+    GetBox(pmin, pmax);
+    auto eps = 1e-8 * (pmax-pmin).Length();
+    for(const auto& se : surfelements)
+      {
+        if(GetBCName(se.index-1) != s1)
+          continue;
 
+        for(const auto& pi : se.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);
+            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);
+          }
+      }
+    return nr;
+  }
 
   void Mesh :: InitPointCurve(double red, double green, double blue) const
   {

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

libsrc/meshing/python_mesh.cpp

@@ -855,6 +855,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::arg("pid2"),
            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)
           {