From aa8dbac6be4c6dce873accf32dcd72774b15b059 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= <joachim.schoeberl@tuwien.ac.at>
Date: Fri, 9 Mar 2018 15:29:50 +0100
Subject: [PATCH] export Trafo to py, IdentifyPeriodic with Trafo

---
 libsrc/csg/identify.cpp        | 19 +++++++++++--------
 libsrc/csg/identify.hpp        |  5 ++++-
 libsrc/csg/python_csg.cpp      |  9 ++++++---
 libsrc/gprim/transform3d.hpp   |  9 ++++++++-
 libsrc/meshing/python_mesh.cpp |  9 ++++++++-
 5 files changed, 37 insertions(+), 14 deletions(-)

diff --git a/libsrc/csg/identify.cpp b/libsrc/csg/identify.cpp
index 2a4da95b..a50c3b5d 100644
--- a/libsrc/csg/identify.cpp
+++ b/libsrc/csg/identify.cpp
@@ -120,9 +120,11 @@ PeriodicIdentification ::
 PeriodicIdentification (int anr,
 			const CSGeometry & ageom,
 			const Surface * as1,
-			const Surface * as2)
-  : Identification(anr, ageom)
+			const Surface * as2,
+                        Transformation<3> atrafo)
+  : Identification(anr, ageom), trafo(atrafo)
 {
+  inv_trafo = trafo.CalcInverse();
   s1 = as1;
   s2 = as2;
 }
@@ -256,25 +258,26 @@ GetIdentifiedPoint (class Mesh & mesh,  int pi)
   const Surface *snew;
   const Point<3> & p = mesh.Point (pi);
 
+  Point<3> hp = p;
   if (s1->PointOnSurface (p))
     {
       snew = s2;
+      hp = trafo(hp);
     }
   else
     {
       if (s2->PointOnSurface (p))
 	{
 	  snew = s1;
+          hp = inv_trafo(hp);
 	}
       else
 	{
-	  cerr << "GetIdenfifiedPoint: Not possible" << endl;
-	  exit (1);
+          throw NgException("GetIdenfifiedPoint: Not possible");
 	}    
     }
   
   // project to other surface
-  Point<3> hp = p;
   snew->Project (hp);
 
   int newpi = 0;
@@ -306,15 +309,15 @@ GetIdentifiedPoint (class Mesh & mesh,  int pi)
 
 void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh)
 {
-  int i, j;
-  for (i = 1; i <= mesh.GetNP(); i++)
+  for (int i = 1; i <= mesh.GetNP(); i++)
     {
       Point<3> p = mesh.Point(i);
       if (s1->PointOnSurface (p))
 	{
 	  Point<3> pp = p;
+          pp = trafo(pp);
 	  s2->Project (pp);
-	  for (j = 1; j <= mesh.GetNP(); j++)
+	  for (int j = 1; j <= mesh.GetNP(); j++)
 	    if (Dist2(mesh.Point(j), pp) < 1e-6)
 	      {
 		mesh.GetIdentifications().Add (i, j, nr);
diff --git a/libsrc/csg/identify.hpp b/libsrc/csg/identify.hpp
index ef3f7556..5211d155 100644
--- a/libsrc/csg/identify.hpp
+++ b/libsrc/csg/identify.hpp
@@ -78,11 +78,14 @@ namespace netgen
   {
     const Surface * s1;
     const Surface * s2;
+    Transformation<3> trafo; // from s1 to s2
+    Transformation<3> inv_trafo; // from s2 to s1
   public:
     PeriodicIdentification (int anr,
 			    const CSGeometry & ageom,
 			    const Surface * as1,
-			    const Surface * as2);
+			    const Surface * as2,
+                            Transformation<3> atrafo = Vec<3>(0,0,0));
     virtual ~PeriodicIdentification ();
     virtual void Print (ostream & ost) const;
     virtual void GetData (ostream & ost) const;
diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp
index d6d60e75..b390fe7f 100644
--- a/libsrc/csg/python_csg.cpp
+++ b/libsrc/csg/python_csg.cpp
@@ -555,7 +555,8 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
          )
     
     .def("PeriodicSurfaces", FunctionPointer
-         ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
+         ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2,
+              Transformation<3> trafo)
           {
             Array<int> si1, si2;
             s1->GetSolid()->GetSurfaceIndices (si1);
@@ -564,9 +565,11 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
             self.AddIdentification 
               (new PeriodicIdentification 
                (self.GetNIdentifications()+1, self, 
-                self.GetSurface (si1[0]), self.GetSurface (si2[0])));
+                self.GetSurface (si1[0]), self.GetSurface (si2[0]),
+                trafo));
           }),
-         py::arg("solid1"), py::arg("solid2")
+         py::arg("solid1"), py::arg("solid2"),
+         py::arg("trafo")=Transformation<3>(Vec<3>(0,0,0))
          )
 
     .def("AddPoint", [] (CSGeometry & self, Point<3> p, int index) -> CSGeometry&
diff --git a/libsrc/gprim/transform3d.hpp b/libsrc/gprim/transform3d.hpp
index e9736ff6..7c966801 100644
--- a/libsrc/gprim/transform3d.hpp
+++ b/libsrc/gprim/transform3d.hpp
@@ -129,7 +129,14 @@ public:
   }
 
   /// 
-  void CalcInverse (Transformation & inv) const;
+  Transformation CalcInverse () const
+  {
+    Transformation inv;
+    // inv.m = Inv(m);
+    ::netgen::CalcInverse (m, inv.m);
+    inv.v = inv.m * (-v);
+    return inv;
+  }
 
   /// this = ta x tb
   void Combine (const Transformation & ta, const Transformation & tb)
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 63f577c9..f4c461ab 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -86,11 +86,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
   m.def ("Pnt", FunctionPointer
            ([](double x, double y) { return Point<2>(x,y); }));
 
-           
+  /*
+    // duplicated functions ????
   m.def ("Pnt", FunctionPointer
            ([](double x, double y, double z) { return Point<3>(x,y,z); }));
   m.def ("Pnt", FunctionPointer
            ([](double x, double y) { return Point<2>(x,y); }));
+  */
 
   py::class_<Vec<2>> (m, "Vec2d")
     .def(py::init<double,double>())
@@ -117,6 +119,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
   m.def ("Vec", FunctionPointer
            ([] (double x, double y) { return Vec<2>(x,y); }));
 
+  py::class_<Transformation<3>> (m, "Trafo")
+    .def(py::init<Vec<3>>())
+    .def("__call__", [] (Transformation<3> trafo, Point<3> p) { return trafo(p); })
+    // .def ("__str__", &ToString<Transformation<3>>)    
+    ;
 
   m.def ("SetTransformation", FunctionPointer
            ([](int dir, double angle)