From 88ec68be6566d2faca538519cd3f6756a446fbcf Mon Sep 17 00:00:00 2001
From: Michael Neunteufel <michael.neunteufel@live.de>
Date: Mon, 24 Jun 2019 19:29:38 +0200
Subject: [PATCH] OptimizeMesh2d works now with curved TRIG6

---
 libsrc/meshing/meshclass.cpp   | 44 ++++++++++++++++++++++++++++++----
 libsrc/meshing/meshfunc2d.cpp  | 17 +++++++++++++
 libsrc/meshing/secondorder.cpp |  8 +++++++
 3 files changed, 65 insertions(+), 4 deletions(-)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index aa077071..1137c9c7 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -3324,8 +3324,8 @@ namespace netgen
     for (int i = 0; i < segments.Size(); i++)
       {
         const Segment & seg = segments[i];
-        pused.Set (seg[0]);
-        pused.Set (seg[1]);
+        for (int j = 0; j < seg.GetNP(); j++)
+          pused.Set (seg[j]);
       }
 
     for (int i = 0; i < openelements.Size(); i++)
@@ -3389,8 +3389,8 @@ namespace netgen
     for (int i = 0; i < segments.Size(); i++)
       {
         Segment & seg = segments[i];
-        seg[0] = op2np[seg[0]];
-        seg[1] = op2np[seg[1]];
+        for (int j = 0; j < seg.GetNP(); j++)
+          seg[j] = op2np[seg[j]];
       }
 
     for (int i = 1; i <= openelements.Size(); i++)
@@ -4789,6 +4789,42 @@ namespace netgen
             //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl;
             //(*testout) << "sol " << sol << endl;
 
+            if (SurfaceElement(element).GetType() ==TRIG6)
+              {
+                netgen::Point<2> lam(1./3,1./3);
+                Vec<3> rhs;
+                Vec<2> deltalam;
+                netgen::Point<3> x;
+                Mat<3,2> Jac,Jact;
+                
+                double delta=1;
+                
+                bool retval;
+                
+                int i = 0;
+                
+                const int maxits = 30;
+                while(delta > 1e-16 && i<maxits)
+                  {
+                    curvedelems->CalcSurfaceTransformation(lam,element-1,x,Jac);
+                    rhs = p-x;
+                    Jac.Solve(rhs,deltalam);
+                    
+                    lam += deltalam;
+                    
+                    delta = deltalam.Length2();
+                    
+                    i++;
+                    //(*testout) << "pcie i " << i << " delta " << delta << " p " << p << " x " << x << " lam " << lam << endl;
+                    //<< "Jac " << Jac << endl;
+                  }
+                
+                if(i==maxits)
+                  return false;
+                
+                sol.X() = lam(0);
+                sol.Y() = lam(1);
+              }
             if (sol.X() >= -eps && sol.Y() >= -eps && 
                 sol.X() + sol.Y() <= 1+eps)
               {
diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp
index 765bbf6a..a5b0afef 100644
--- a/libsrc/meshing/meshfunc2d.cpp
+++ b/libsrc/meshing/meshfunc2d.cpp
@@ -10,6 +10,16 @@ namespace netgen
 
     mesh.CalcSurfacesOfNode();
 
+    bool secondorder = mesh.GetNP() > mesh.GetNV();
+
+
+    if (secondorder)
+      {
+      for (SurfaceElementIndex ei = 0; ei < mesh.GetNSE(); ei++)
+        mesh[ei].SetType(TRIG);
+      }
+    mesh.Compress();
+
     const char * optstr = mp.optimize2d.c_str();
     int optsteps = mp.optsteps2d;
 
@@ -51,6 +61,13 @@ namespace netgen
 	      cerr << "Optimization code " << optstr[j-1] << " not defined" << endl;
 	    }  
 	}
+    if (secondorder)
+      {
+        if (mesh.GetGeometry())
+          mesh.GetGeometry()->GetRefinement().MakeSecondOrder(mesh);
+        else
+          Refinement().MakeSecondOrder(mesh);
+      }
   }
 
 }
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 49f91662..c70841b3 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -23,6 +23,14 @@ namespace netgen
 
     INDEX_2_HASHTABLE<PointIndex> between(mesh.GetNP() + 5);
 
+    for (SegmentIndex si = 0; si < mesh.GetNSeg(); si++)
+      {
+        auto & seg = mesh[si];
+        if (seg.GetType() == SEGMENT3)
+          between.Set(INDEX_2::Sort(seg[0],seg[1]), seg[2]);
+      }
+    
+    
     for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
       {
 	const Element2d & el = mesh[sei];