FindElementOfPoint<1> for 2d meshes for curved segments

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];
@@ -1059,6 +1067,7 @@ namespace netgen
                 }
             }
             }
+        }
         break;
       case 3:
       default:

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<int> edges;
+        topology.GetSurfaceElementEdges(velement, edges);
+        Array<SegmentIndex> 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

libsrc/meshing/topology.cpp

@@ -620,7 +620,8 @@ namespace netgen
             ned += hv;
           }
         edge2vert.SetSize(ned);
-
+        edge2segment.SetSize(ned);
+        edge2segment = -1;
 
         // INDEX_CLOSED_HASHTABLE<int> v2eht(2*max_edge_on_vertex+10);
 	// NgArray<int> vertex2;
@@ -689,6 +690,7 @@ namespace netgen
                                       break;
                                     case 1:
                                       segedges[elnr].nr = edgenum;
+                                      edge2segment[edgenum] = elnr;
                                       // segedges[elnr].orient = edgedir;
                                       break;
                                     }

libsrc/meshing/topology.hpp

@@ -63,6 +63,7 @@ class MeshTopology
   NgArray<T_FACE> surffaces;
   NgArray<INDEX_2> surf2volelement;
   NgArray<int> face2surfel;
+  Array<SegmentIndex> edge2segment;
   TABLE<ElementIndex,PointIndex::BASE> vert2element;
   TABLE<SurfaceElementIndex,PointIndex::BASE> vert2surfelement;
   TABLE<SegmentIndex,PointIndex::BASE> vert2segment;
@@ -170,6 +171,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<ElementIndex> & elements) const;
   NgFlatArray<ElementIndex> GetVertexElements (int vnr) const
   { return vert2element[vnr]; }