diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index bb93df5d..072bf8f0 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1089,7 +1089,7 @@ namespace netgen if (ind > 0) { - if(mesh->SurfaceElement(ind).GetType()==QUAD) + if(mesh->SurfaceElement(ind).GetType()==QUAD || mesh->SurfaceElement(ind).GetType()==TRIG6) { lami[0] = lam3[0]; lami[1] = lam3[1]; diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index aa077071..f5dfc7de 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -123,6 +123,8 @@ namespace netgen if (mesh2.cd3names[i]) cd3names[i] = new string(*mesh2.cd3names[i]); else cd3names[i] = 0; + numvertices = mesh2.numvertices; + return *this; } @@ -3324,8 +3326,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 +3391,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 +4791,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 && iCalcSurfaceTransformation(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 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];