diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index e417be3c..10e8c498 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -1035,6 +1035,14 @@ namespace netgen case 2: { Point<3> p(hp[0], hp[1],0); + try + { + auto ind = mesh->GetSurfaceElementOfPoint(p, lami, nullptr, + build_searchtree); + return ind - 1; + } + catch(NgException e) // quads not implemented curved yet + { for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) { auto & seg = (*mesh)[si]; @@ -1058,6 +1066,7 @@ namespace netgen return si; } } + } } break; case 3: diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 9ca15f76..4f18c139 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5268,7 +5268,7 @@ namespace netgen //(*testout) << "col1 " << col1 << " col2 " << col2 << " col3 " << col3 << " rhs " << rhs << endl; //(*testout) << "sol " << sol << endl; - if (SurfaceElement(element).GetType() ==TRIG6) + if (SurfaceElement(element).GetType() ==TRIG6 || curvedelems->IsSurfaceElementCurved(element-1)) { netgen::Point<2> lam(1./3,1./3); Vec<3> rhs; @@ -5679,6 +5679,47 @@ namespace netgen { if (dimension == 2) { + double vlam[3]; + int velement = GetElementOfPoint(p, vlam, NULL, build_searchtree, allowindex); + if(velement == 0) + return 0; + + vlam[2] = 1.-vlam[0] - vlam[1]; + NgArray edges; + topology.GetSurfaceElementEdges(velement, edges); + Array segs(edges.Size()); + for(auto i : Range(edges)) + segs[i] = topology.GetSegmentOfEdge(edges[i]); + + for(auto i : Range(segs)) + { + if(IsInvalid(segs[i])) + continue; + auto& el = SurfaceElement(velement); + if(el.GetType() == TRIG) + { + double seg_lam; + double lam; + auto seg = LineSegment(segs[i]); + for(auto k : Range(3)) + { + if(seg[0] == el[k]) + lam = vlam[k]; + if(seg[1] == el[k]) + seg_lam = vlam[k]; + } + if(1.- seg_lam - lam < 1e-5) + { + // found point close to segment -> use barycentric coordinates directly + lami[0] = lam; + return int(segs[i])+1; + } + } + else + throw NgException("Quad not implemented yet!"); + } + + return 0; throw NgException("GetSurfaceElementOfPoint not yet implemented for 2D meshes"); } else diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index ac571be7..0174d892 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -620,7 +620,8 @@ namespace netgen ned += hv; } edge2vert.SetSize(ned); - + edge2segment.SetSize(ned); + edge2segment = -1; // INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); // NgArray vertex2; @@ -689,6 +690,7 @@ namespace netgen break; case 1: segedges[elnr].nr = edgenum; + edge2segment[edgenum] = elnr; // segedges[elnr].orient = edgedir; break; } diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 57205600..53fbc0ff 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -63,6 +63,7 @@ class MeshTopology NgArray surffaces; NgArray surf2volelement; NgArray face2surfel; + Array edge2segment; TABLE vert2element; TABLE vert2surfelement; TABLE vert2segment; @@ -169,6 +170,8 @@ public: } int GetFace2SurfaceElement (int fnr) const { return face2surfel[fnr-1]; } + + SegmentIndex GetSegmentOfEdge(int edgenr) const { return edge2segment[edgenr-1]; } void GetVertexElements (int vnr, NgArray & elements) const; NgFlatArray GetVertexElements (int vnr) const