From 2f18c2b1f72660174ab5d7643d1babc8d5826533 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Tue, 1 Sep 2020 20:50:03 +0200
Subject: [PATCH] Mesh::Mirror()

---
 libsrc/meshing/meshclass.cpp   | 86 ++++++++++++++++++++++++++++++++++
 libsrc/meshing/meshclass.hpp   |  1 +
 libsrc/meshing/python_mesh.cpp |  1 +
 3 files changed, 88 insertions(+)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 6bcd6e0c..d84db555 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6860,4 +6860,90 @@ namespace netgen
     if (surfelementht)
       surfelementht->PrintMemInfo (cout);
   }
+
+  shared_ptr<Mesh> Mesh :: Mirror ( netgen::Point<3> p_plane, Vec<3> n_plane )
+  {
+    Mesh & m = *this;
+    auto nm_ = make_shared<Mesh>();
+    Mesh & nm = *nm_;
+    nm = m;
+
+    Point3d pmin, pmax;
+    GetBox(pmin, pmax);
+    auto v = pmax-pmin;
+    double eps = v.Length()*1e-8;
+
+    auto onPlane = [&] (const MeshPoint & p) -> bool
+    {
+      auto v = p_plane-p;
+      auto l = v.Length();
+      if(l<eps) return true;
+
+      auto ret = fabs(v*n_plane)/l;
+      return fabs(v*n_plane) < eps;
+    };
+
+    auto mirror = [&] (PointIndex pi) -> PointIndex
+    {
+      auto & p = m[pi];
+
+      auto v = p_plane-p;
+      auto l = v.Length();
+      if(l<eps)
+        return pi;
+
+      if(fabs(v*n_plane)/l < eps)
+        return pi;
+
+      auto new_point = p + 2*(v*n_plane)*n_plane;
+      return nm.AddPoint( new_point, p.GetLayer(), p.Type() );
+    };
+
+    Array<PointIndex, PointIndex> point_map;
+    point_map.SetSize(GetNP());
+    point_map = -1;
+
+    for(auto pi : Range(points))
+      point_map[pi] = mirror(pi);
+
+    for(auto & el : VolumeElements())
+    {
+      auto nel = el;
+      for(auto i : Range(el.GetNP()))
+        nel[i] = point_map[el[i]];
+      nm.AddVolumeElement(nel);
+    }
+
+    for (auto ei : Range(SurfaceElements()))
+    {
+      auto & el = m[ei];
+      auto nel = el;
+      for(auto i : Range(el.GetNP()))
+        nel[i] = point_map[el[i]];
+
+      if(!(nel==el))
+        nm.AddSurfaceElement(nel);
+    }
+
+    for (auto ei : Range(LineSegments()))
+    {
+      auto & el = LineSegments()[ei];
+      auto nel = el;
+      bool is_same = true;
+
+      for(auto i : Range(el.GetNP()))
+      {
+        auto pi = el[i];
+        nel[i] = point_map[pi];
+        if(point_map[pi]!=pi)
+          is_same = false;
+      }
+
+      if(!is_same)
+        nm.AddSegment(nel);
+    }
+
+    return nm_;
+  }
+
 }
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 5139cb8c..59ed7b85 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -923,6 +923,7 @@ namespace netgen
       NgArray<int> & segment_weights){ }
 #endif
 
+    shared_ptr<Mesh> Mirror( netgen::Point<3> p, Vec<3> n );
 
   };
 
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 07856c1c..5b2ea5c1 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1150,6 +1150,7 @@ grow_edges : bool = False
           })
     .def ("CalcTotalBadness", &Mesh::CalcTotalBad)
     .def ("GetQualityHistogram", &Mesh::GetQualityHistogram)
+    .def("Mirror", &Mesh::Mirror);
     ;
 
   m.def("ImportMesh", [](const string& filename)