Merge branch 'fix_getelementofpoint' into 'master'

Restructure code in GetElementOfPoint

See merge request jschoeberl/netgen!424

libsrc/meshing/meshclass.cpp

@@ -5609,52 +5609,8 @@ namespace netgen
       return -1;
 
     if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE()))
-      {
+      return GetSurfaceElementOfPoint(p, lami, indices, build_searchtree, allowindex);
-        int ne;
-        int ps_startelement = 0;  // disable global buffering
 
-        if(ps_startelement != 0 && ps_startelement <= GetNSE() && PointContainedIn2DElement(p,lami,ps_startelement))
-          return ps_startelement;
-
-        NgArray<int> locels;
-        if (elementsearchtree || build_searchtree)
-          {
-            // update if necessary:
-            const_cast<Mesh&>(*this).BuildElementSearchTree ();
-            // double tol = elementsearchtree->Tolerance();
-            // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol);
-            // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol);
-            elementsearchtree->GetIntersecting (p, p, locels);
-            ne = locels.Size();
-          }
-        else
-          ne = GetNSE();
-
-        for (int i = 1; i <= ne; i++)
-          {
-            int ii;
-
-            if (elementsearchtree)
-              ii = locels.Get(i);
-            else
-              ii = i;
-
-            if(ii == ps_startelement) continue;
-
-            if(indices != NULL && indices->Size() > 0)
-              {
-                bool contained = indices->Contains(SurfaceElement(ii).GetIndex());
-                if((allowindex && !contained) || (!allowindex && contained)) continue;
-              }
-
-            if(PointContainedIn2DElement(p,lami,ii)) return ii;
-
-          }
-        return 0;
-      }
-    else
-
-      {
     int ps_startelement = 0;  // disable global buffering
     // int i, j;
     int ne;
@@ -5723,11 +5679,8 @@ namespace netgen
             return ii;
           }
       }
-
-
     return 0;
   }
-  }
 
 
 
@@ -5809,14 +5762,9 @@ namespace netgen
         //(*testout) << "p " << p << endl;
         //(*testout) << "velement " << velement << endl;
 
-        if (!GetNE() && GetNSE() )
+        // first try to find a volume element containing p and project to face
+        if(velement!=0)
         {
-            lami[0] = vlam[0];
-            lami[1] = vlam[1];
-            lami[2] = vlam[2];
-            return velement;
-          }
-        
           NgArray<int> faces;
           topology.GetElementFaces(velement,faces);
 
@@ -5831,6 +5779,9 @@ namespace netgen
             {
               if(faces[i] == 0)
                 continue;
+              auto sel = SurfaceElement(faces[i]);
+              if(indices && indices->Size() != 0 && !indices->Contains(sel.GetIndex()))
+                continue;
 
               auto & el = VolumeElement(velement);
 
@@ -5842,8 +5793,6 @@ namespace netgen
                 {
                   // found volume point very close to a face -> use barycentric coordinates directly
                   lami[2] = 0.0;
-                auto sel = SurfaceElement(faces[i]);
-
                   for(auto j : Range(1,3))
                     for(auto k : Range(4))
                       if(sel[j] == el[k])
@@ -5852,33 +5801,47 @@ namespace netgen
                 }
               }
 
-            if(indices && indices->Size() != 0)
-              {
-                if(indices->Contains(SurfaceElement(faces[i]).GetIndex()) &&
-                   PointContainedIn2DElement(p,lami,faces[i],true))
-                  return faces[i];
-              }
-            else
-              {
               if(PointContainedIn2DElement(p,lami,faces[i],true))
-                  {
-                    //(*testout) << "found point " << p << " in sel " << faces[i]
-                    //	       << ", lam " << lami[0] << ", " << lami[1] << ", " << lami[2] << endl;
                 return faces[i];
             }
         }
-          }
 
-        NgArray<int> faces2;
+        // Did't find any matching face of a volume element, search 2d elements directly
-        topology.GetElementFaces(velement,faces2);
+        int ne;
-        /*
+
-        cout << "no matching surf element" << endl
+        NgArray<int> locels;
-             << "p = " << p << endl
+        // TODO: build search tree for surface elements
-             << "faces-orig = " << faces2 << endl
+        if (!GetNE() && GetNSE() && (elementsearchtree || build_searchtree))
-             << "faces = " << faces << endl
+          {
-             << ", vol el = " << velement
+            // update if necessary:
-             << ", vlam = " << vlam[0] << "," << vlam[1] << "," << vlam[2] << endl;
+            const_cast<Mesh&>(*this).BuildElementSearchTree ();
-        */
+            // double tol = elementsearchtree->Tolerance();
+            // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol);
+            // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol);
+            elementsearchtree->GetIntersecting (p, p, locels);
+            ne = locels.Size();
+          }
+        else
+          ne = GetNSE();
+
+        for (int i = 1; i <= ne; i++)
+          {
+            int ii;
+
+            if (locels.Size())
+              ii = locels.Get(i);
+            else
+              ii = i;
+
+            if(indices != NULL && indices->Size() > 0)
+              {
+                bool contained = indices->Contains(SurfaceElement(ii).GetIndex());
+                if((allowindex && !contained) || (!allowindex && contained)) continue;
+              }
+
+            if(PointContainedIn2DElement(p,lami,ii)) return ii;
+
+          }
       }
 
     return 0;