From 95301e11ba966eb50c9deb7648b74cd1def7255a Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Fri, 15 Oct 2021 18:52:11 +0200
Subject: [PATCH] mesh.SecondOrder : fix Segment mapping

---
 libsrc/gprim/geomobjects.hpp   |  2 +-
 libsrc/gprim/geomops.hpp       |  9 +++++++++
 libsrc/meshing/curvedelems.cpp | 25 +++++++++++++++++++++++++
 libsrc/meshing/secondorder.cpp |  1 +
 4 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/libsrc/gprim/geomobjects.hpp b/libsrc/gprim/geomobjects.hpp
index f0f8edb4..5cfc2eb9 100644
--- a/libsrc/gprim/geomobjects.hpp
+++ b/libsrc/gprim/geomobjects.hpp
@@ -43,7 +43,7 @@ namespace netgen
     Point (const Point<D,T2> & p2)
     { for (int i = 0; i < D; i++) x[i] = p2(i); }
 
-    explicit Point (const Vec<D> & v)
+    explicit Point (const Vec<D,T> & v)
     { for (int i = 0; i < D; i++) x[i] = v(i); }
 
 
diff --git a/libsrc/gprim/geomops.hpp b/libsrc/gprim/geomops.hpp
index 6b06cb6e..fcb127ed 100644
--- a/libsrc/gprim/geomops.hpp
+++ b/libsrc/gprim/geomops.hpp
@@ -78,6 +78,15 @@ namespace netgen
     return res;
   }
 
+  template <int D, int SW>
+  inline Vec<D, SIMD<double,SW>> operator* (SIMD<double,SW> s, Vec<D,double> b)
+  {
+    Vec<D,SIMD<double,SW>> res;
+    for (int i = 0; i < D; i++)
+      res(i) = s * b(i);
+    return res;
+  }
+  
 
   template <int D>
   inline double operator* (Vec<D> a, Vec<D> b)
diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp
index 79f031a4..21237782 100644
--- a/libsrc/meshing/curvedelems.cpp
+++ b/libsrc/meshing/curvedelems.cpp
@@ -1429,6 +1429,31 @@ namespace netgen
     info.order = order;
     info.ndof = info.nv = 2;
 
+    if (order == 1)
+      {
+        auto & seg = mesh.LineSegment(elnr);
+        if (seg.GetNP() == 2)
+          {
+            if (x)
+              *x = Point<3,T>((1-xi) * Vec<3>(mesh[seg.PNums()[0]]) + xi * Vec<3>(mesh[seg.PNums()[1]]));
+            if (dxdxi)
+              *dxdxi = Vec<3>(mesh[seg.PNums()[1]])-Vec<3>(mesh[seg.PNums()[0]]);
+          }
+        else
+          {
+            if (x)
+              {
+                *x = Point<3,T>(2*(xi-1)*(xi-0.5) * Vec<3>(mesh[seg.PNums()[1]])
+                                + 4*xi*(1-xi) * Vec<3>(mesh[seg.PNums()[2]])
+                                + 2*xi*(xi-0.5) * Vec<3>(mesh[seg.PNums()[0]]));
+              }
+            if (dxdxi)
+              *dxdxi = T(1.0) * (Vec<3>(mesh[seg.PNums()[1]])-Vec<3>(mesh[seg.PNums()[0]]))
+                + (4-8*xi)*Vec<3>(mesh[seg.PNums()[2]]);
+          }
+        return;
+      }
+    
     if (info.order > 1)
       {
 	const MeshTopology & top = mesh.GetTopology();
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 3a7368a9..9e03d5fc 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -110,6 +110,7 @@ namespace netgen
 				   EDGEPOINT);
 	    between.Set (i2, el[2]);
 	  }
+        el.SetCurved(true);
       }
 
     // refine surface elements