Merge branch 'findpointinelement' into 'master'

Improvements to FindPointInElement interface code

See merge request ngsolve/netgen!701

libsrc/interface/nginterface_v2.cpp

@@ -1011,68 +1011,19 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    switch (mesh->GetDimension())
+    Point<3> p(hp[0], 0,0);
+    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
       {
-      case 1:
+        auto & seg = (*mesh)[si];
-        {
+        Point<3> p1 = (*mesh)[seg[0]];
-          Point<3> p(hp[0], 0,0);
+        Point<3> p2 = (*mesh)[seg[1]];
-          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
+        double lam = (p(0)-p1(0)) / (p2(0)-p1(0));
-            {
+        if (lam >= -1e-10 && lam <= 1+1e-10)
-              auto & seg = (*mesh)[si];
+          {
-              Point<3> p1 = (*mesh)[seg[0]];
+            lami[0] = 1-lam;
-              Point<3> p2 = (*mesh)[seg[1]];
+            return si;
-              double lam = (p(0)-p1(0)) / (p2(0)-p1(0));
+          }
-              if (lam >= -1e-10 && lam <= 1+1e-10)
-                {
-                  lami[0] = 1-lam;
-                  return si;
-                }
-            }
-        }
-        break;
-      case 2:
-        {
-          Point<3> p(hp[0], hp[1],0);
-          try
-            {
-              auto ind = mesh->GetSurfaceElementOfPoint(p, lami, nullptr,
-                                                        build_searchtree);
-              return ind - 1;
-            }
-          catch(const NgException & e) // quads not implemented curved yet
-            {
-          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
-            {
-              auto & seg = (*mesh)[si];
-              Point<3> p1 = (*mesh)[seg[0]];
-              Point<3> p2 = (*mesh)[seg[1]];
-              double lam;
-              double r;
-              if (fabs(p2[0]-p1[0]) >= fabs(p2[1]-p1[1]))
-                {
-                  lam = (p[0]-p1[0])/(p2[0]-p1[0]);
-                  r = p[1] - p1[1] - lam*(p2[1]-p1[1]);
-                }
-              else
-                {
-                  lam = (p[1]-p1[1])/(p2[1]-p1[1]);
-                  r = p[0] - p1[0] - lam*(p2[0]-p1[0]);
-                }
-              if ( lam >= -1e-10 && lam <= 1+1e-10 && fabs(r) <= 1e-10 )
-                {
-                  lami[0] = 1-lam;
-                  return si;
-                }
-            }
-            }
-        }
-        break;
-      case 3:
-      default:
-        throw Exception("FindElementOfPoint<1> only implemented for mesh-dimension 1 and 2!");
-        break;
       }
- 
     return -1;
   }
 
@@ -1083,37 +1034,13 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    NgArray<int> dummy(numind);
+    if(build_searchtree)
-    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
+      mesh->BuildElementSearchTree(2);
-    
+    Point<3> pp(p[0], p[1], 0.);
-    double lam3[3];
+    if(mesh->GetDimension() == 3)
-    int ind;
+      pp[2] = p[2];
-
+    NgArray<int> ind(numind, indices);
-    if (mesh->GetDimension() == 2)
+    return Find2dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(2));
-      {
-        Point<3> p2d(p[0], p[1], 0);
-        ind = mesh->GetElementOfPoint(p2d, lam3, &dummy, build_searchtree);
-      }
-    else
-      {
-        Point3d p3d(p[0], p[1], p[2]);
-        ind = mesh->GetSurfaceElementOfPoint(p3d, lam3, &dummy, build_searchtree);
-      }
-    
-    if (ind > 0)
-      {
-        if(mesh->SurfaceElement(ind).GetType()==QUAD || mesh->SurfaceElement(ind).GetType()==TRIG6)
-          {
-            lami[0] = lam3[0];
-            lami[1] = lam3[1];
-          }
-        else 
-          {
-            lami[0] = 1-lam3[0]-lam3[1];
-            lami[1] = lam3[0];
-          }
-      }
-    return ind-1;
   }
 
 
@@ -1124,13 +1051,11 @@ namespace netgen
    int * const indices, int numind) const
 
   {
-    NgArray<int> dummy(numind);
+    if(build_searchtree)
-    for (int i = 0; i < numind; i++) dummy[i] = indices[i]+1;
+      mesh->BuildElementSearchTree(3);
-    
+    Point<3> pp(p[0], p[1], p[2]);
-    Point<3> p3d(p[0], p[1], p[2]);
+    NgArray<int> ind(numind, indices);
-    int ind = 
+    return Find3dElement(*mesh, pp, lami, &ind, mesh->GetElementSearchTree(3));
-      mesh->GetElementOfPoint(p3d, lami, &dummy, build_searchtree);
-    return ind-1;
   }
 
   void Ngx_Mesh :: Curve (int order)

libsrc/interface/writeOpenFOAM15x.cpp

@@ -607,7 +607,7 @@ namespace netgen
       const_cast<Mesh&> (mesh).Compress();
       const_cast<Mesh&> (mesh).CalcSurfacesOfNode();
       const_cast<Mesh&> (mesh).RebuildSurfaceElementLists();
-      const_cast<Mesh&> (mesh).BuildElementSearchTree();
+      const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
 
 
       int np = mesh.GetNP();

libsrc/interface/writefluent.cpp

@@ -79,7 +79,7 @@ void WriteFluentFormat (const Mesh & mesh,
   snprintf (str, size(str), "(13 (4 1 %x 2 3)(",noverbface); //hexadecimal!!!
   outfile << str << endl;
 
-  const_cast<Mesh&> (mesh).BuildElementSearchTree();
+  const_cast<Mesh&> (mesh).BuildElementSearchTree(3);
 
   for (i = 1; i <= ne; i++)
     {

libsrc/meshing/meshclass.cpp

@@ -17,7 +17,7 @@ namespace netgen
 		     double * lami,
 		     const NgArray<int> * const indices,
 		     BoxTree<3> * searchtree,
-		     const bool allowindex = true)
+		     const bool allowindex)
   {
     int ne = 0;
     NgArray<int> locels;
@@ -79,13 +79,13 @@ namespace netgen
 		     double * lami,
 		     const NgArray<int> * const indices,
 		     BoxTree<3> * searchtree,
-		     const bool allowindex = true)
+		     const bool allowindex)
   {
     double vlam[3];
     int velement = 0;
 
     if(mesh.GetNE())
-      velement = Find3dElement(mesh, p,vlam,NULL,searchtree,allowindex);
+      velement = Find3dElement(mesh, p,vlam,NULL,searchtree ? mesh.GetElementSearchTree(3) : nullptr,allowindex);
 
     //(*testout) << "p " << p << endl;
     //(*testout) << "velement " << velement << endl;
@@ -141,7 +141,7 @@ namespace netgen
 
     NgArray<int> locels;
     // TODO: build search tree for surface elements
-    if (!mesh.GetNE() && searchtree)
+    if (searchtree)
       {
         searchtree->GetIntersecting (p, p, locels);
         ne = locels.Size();
@@ -235,8 +235,8 @@ namespace netgen
     : topology(*this), surfarea(*this)
   {
     lochfunc = {nullptr};
-    elementsearchtree = nullptr;
+    for(auto i : Range(4))
-    elementsearchtreets = NextTimeStamp();
+      elementsearchtreets[i] = NextTimeStamp();
     majortimestamp = timestamp = NextTimeStamp();
     hglob = 1e10;
     hmin = 0;
@@ -5251,122 +5251,91 @@ namespace netgen
     RebuildSurfaceElementLists();
   }
 
-  void Mesh :: BuildElementSearchTree ()
+  void Mesh :: BuildElementSearchTree (int dim)
   {
-    if (elementsearchtreets == GetTimeStamp()) return;
+    if (elementsearchtreets[dim] == GetTimeStamp())
+      return;
 
     {
       std::lock_guard<std::mutex> guard(buildsearchtree_mutex);
-      if (elementsearchtreets != GetTimeStamp())
+      // check again to see if some other thread built while waiting for lock
+      if (elementsearchtreets[dim] == GetTimeStamp()) return;
+          
+      PrintMessage (4, "Rebuild element searchtree dim " + ToString(dim));
+          
+
+      Point3d pmin, pmax;
+      GetBox(pmin, pmax);
+      Box<3> box(pmin, pmax);
+      if (dim > 1)
+        elementsearchtree[dim] = make_unique<BoxTree<3>>(box);
+      else
+        elementsearchtree[dim] = nullptr; // not yet implemented
+
+      if (dim == 3)
         {
-          NgLock lock(mutex);
+          for(auto ei : volelements.Range())
-          lock.Lock();
-          
-          PrintMessage (4, "Rebuild element searchtree");
-          
-          elementsearchtree = nullptr;
-          
-          int ne = (dimension == 2) ? GetNSE() : GetNE();
-          if (dimension == 3 && !GetNE() && GetNSE())
-            ne = GetNSE();
-
-          if (ne) 
             {
-              if (dimension == 2 || (dimension == 3 && !GetNE()) )
+              const auto& el = volelements[ei];
+              Box<3> box (Box<3>::EMPTY_BOX);
+              for (auto pi : el.PNums())
+                box.Add (points[pi]);
+
+              if(el.IsCurved() && curvedelems->IsElementCurved(ei))
                 {
-                  Box<3> box (Box<3>::EMPTY_BOX);
+                  // add edge/face midpoints to box
-                  for (SurfaceElementIndex sei = 0; sei < ne; sei++)
+                  auto eltype = el.GetType();
-                    // box.Add (points[surfelements[sei].PNums()]);
+                  const auto verts = topology.GetVertices(eltype);
-                    for (auto pi : surfelements[sei].PNums())
-                      box.Add (points[pi]);
-                  
-                  box.Increase (1.01 * box.Diam());
-                  elementsearchtree = make_unique<BoxTree<3>> (box);
-                  
-                  for (SurfaceElementIndex sei = 0; sei < ne; sei++)
-                    {
-                      //  box.Set (points[surfelements[sei].PNums()]);
 
-                      Box<3> box (Box<3>::EMPTY_BOX);
+                  const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
-                      for (auto pi : surfelements[sei].PNums())
+                  for (const auto & edge: edges) {
-                        box.Add (points[pi]);
+                    netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
+                    auto p = netgen::Point<3>(0.0);
+                    curvedelems->CalcElementTransformation(lam,ei,p);
+                    box.Add(p);
+                  }
 
-                      auto & el = surfelements[sei];
+                  const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
-                      if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(sei))
+                  for (const auto & face: faces) {
-                        {
+                    netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
-                          netgen::Point<2>  lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
+                    if(face[3] != -1) {
-                          for (auto lam : lami)
+                      lam += netgen::Vec<3>(verts[face[3]]);
-                            {
+                      lam *= 0.25;
-                              netgen::Point<3> x;
-                              Mat<3,2> Jac;
-                          
-                              curvedelems->CalcSurfaceTransformation(lam,sei,x,Jac);
-                              box.Add (x);
-                            }
-                          box.Scale(1.2);
-                        }
-                      elementsearchtree -> Insert (box, sei+1);
                     }
+                    else
+                      lam *= 1.0/3;
+                    auto p = netgen::Point<3>(0.0);
+                    curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
+                    box.Add(p);
+                  }
                 }
-              else
+              box.Scale(1.2);
+              elementsearchtree[dim] -> Insert (box, ei+1);
+            }
+        }
+      else if (dim == 2)
+        {
+          for (auto ei : Range(surfelements))
+            {
+              const auto& el = surfelements[ei];
+              Box<3> box (Box<3>::EMPTY_BOX);
+              for (auto pi : el.PNums())
+                box.Add (points[pi]);
+
+              if(el.IsCurved() && curvedelems->IsSurfaceElementCurved(ei))
                 {
-                  Box<3> box (Box<3>::EMPTY_BOX);
+                  netgen::Point<2>  lami [4] = {netgen::Point<2>(0.5,0), netgen::Point<2>(0,0.5), netgen::Point<2>(0.5,0.5), netgen::Point<2>(1./3,1./3)};
-                  for (ElementIndex ei = 0; ei < ne; ei++)
+                  for (auto lam : lami)
-                    // box.Add (points[volelements[ei].PNums()]);
-                    for (auto pi : volelements[ei].PNums())
-                      box.Add (points[pi]);
-                  
-                  box.Increase (1.01 * box.Diam());
-                  elementsearchtree = make_unique<BoxTree<3>> (box);
-                  
-                  for (ElementIndex ei = 0; ei < ne; ei++)
                     {
-                      // box.Set (points[volelements[ei].PNums()]);
+                      netgen::Point<3> x;
+                      Mat<3,2> Jac;
 
-                      Box<3> box (Box<3>::EMPTY_BOX);
+                      curvedelems->CalcSurfaceTransformation(lam,ei,x,Jac);
-                      for (auto pi : volelements[ei].PNums())
+                      box.Add (x);
-                        box.Add (points[pi]);
-
-                      auto & el = volelements[ei];
-                      if(el.IsCurved() && curvedelems->IsElementCurved(ei))
-                      {
-                        // add edge/face midpoints to box
-                        auto eltype = el.GetType();
-                        const auto verts = topology.GetVertices(eltype);
-
-
-                        const auto edges = FlatArray<const ELEMENT_EDGE>(topology.GetNEdges(eltype), topology.GetEdges0(eltype));
-                        for (const auto & edge: edges) {
-                          netgen::Point<3> lam = netgen::Point<3>(0.5* (verts[edge[0]] + verts[edge[1]]));
-                          auto p = netgen::Point<3>(0.0);
-                          curvedelems->CalcElementTransformation(lam,ei,p);
-                          box.Add(p);
-                        }
-
-                        const auto faces = FlatArray<const ELEMENT_FACE>(topology.GetNFaces(eltype), topology.GetFaces0(eltype));
-                        for (const auto & face: faces) {
-                          netgen::Vec<3> lam = netgen::Vec<3>(verts[face[0]] + verts[face[1]] + verts[face[2]]);
-                          if(face[3] != -1) {
-                            lam += netgen::Vec<3>(verts[face[3]]);
-                            lam *= 0.25;
-                          }
-                          else
-                            lam *= 1.0/3;
-                          auto p = netgen::Point<3>(0.0);
-                          curvedelems->CalcElementTransformation(netgen::Point<3>(lam),ei,p);
-                          box.Add(p);
-                        }
-
-
-                        box.Scale(1.2);
-                      }
-
-
-                      elementsearchtree -> Insert (box, ei+1);
                     }
+                  box.Scale(1.2);
                 }
-              
+              elementsearchtree[dim] -> Insert (box, ei+1);
-              elementsearchtreets = GetTimeStamp();
             }
         }
     }
@@ -6073,13 +6042,17 @@ namespace netgen
     	 (dimension == 3 && !GetNE() && !GetNSE()) )
       return -1;
 
-    if (build_searchtree)
-      const_cast<Mesh&>(*this).BuildElementSearchTree ();
 
     if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE()))
-      return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
+      {
+        if (build_searchtree)
+          const_cast<Mesh&>(*this).BuildElementSearchTree (2);
+        return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
+      }
 
-    return Find3dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
+    if (build_searchtree)
+      const_cast<Mesh&>(*this).BuildElementSearchTree (3);
+    return Find3dElement(*this, p, lami, indices, elementsearchtree[3].get(), allowindex);
   }
 
 
@@ -6109,13 +6082,14 @@ namespace netgen
                                         bool build_searchtree,
                                         const bool allowindex) const
   {
-    if (!GetNE() && build_searchtree)
-      const_cast<Mesh&>(*this).BuildElementSearchTree ();
-
     if (dimension == 2)
-        return Find1dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
+      return Find1dElement(*this, p, lami, indices, nullptr, allowindex);
     else
-        return Find2dElement(*this, p, lami, indices, elementsearchtree.get(), allowindex);
+      {
+        if (build_searchtree)
+          const_cast<Mesh&>(*this).BuildElementSearchTree(2);
+        return Find2dElement(*this, p, lami, indices, elementsearchtree[2].get(), allowindex);
+      }
     return 0;
   }
 
@@ -6123,7 +6097,7 @@ namespace netgen
   void Mesh::GetIntersectingVolEls(const Point3d& p1, const Point3d& p2, 
                                    NgArray<int> & locels) const
   {
-    elementsearchtree->GetIntersecting (p1, p2, locels);
+    elementsearchtree[3]->GetIntersecting (p1, p2, locels);
   }
 
   void Mesh :: SplitIntoParts()

libsrc/meshing/meshclass.hpp

@@ -148,9 +148,9 @@ namespace netgen
     int numvertices;
 
     /// geometric search tree for interval intersection search
-    unique_ptr<BoxTree<3>> elementsearchtree;
+    unique_ptr<BoxTree<3>> elementsearchtree[4] = {nullptr,nullptr,nullptr,nullptr};
     /// time stamp for tree
-    mutable int elementsearchtreets;
+    mutable size_t elementsearchtreets[4];
 
     /// element -> face, element -> edge etc ...
     MeshTopology topology;
@@ -663,7 +663,11 @@ namespace netgen
 
 
     /// build box-search tree
-    DLL_HEADER void BuildElementSearchTree ();
+    DLL_HEADER void BuildElementSearchTree (int dim);
+    BoxTree<3>* GetElementSearchTree (int dim) const
+    {
+      return elementsearchtree[dim].get();
+    }
 
     void SetPointSearchStartElement(const int el) const {ps_startelement = el;}
 
@@ -1036,6 +1040,20 @@ namespace netgen
   }
 
   DLL_HEADER void AddFacesBetweenDomains(Mesh & mesh);
+
+  int Find2dElement (const Mesh& mesh,
+                     const netgen::Point<3> & p,
+		     double * lami,
+		     const NgArray<int> * const indices,
+		     BoxTree<3> * searchtree,
+		     const bool allowindex = true);
+
+  int Find3dElement (const Mesh& mesh,
+                     const netgen::Point<3> & p,
+		     double * lami,
+		     const NgArray<int> * const indices,
+		     BoxTree<3> * searchtree,
+		     const bool allowindex = true);
   
 }
 

libsrc/meshing/parallelmesh.cpp

@@ -966,7 +966,7 @@ namespace netgen
     self.ident = make_unique<Identifications> (self);
     self.topology = MeshTopology(*this);
     self.topology.Update();
-    self.BuildElementSearchTree();
+    self.BuildElementSearchTree(3);
     
     // const_cast<Mesh&>(*this).DeleteMesh();
 

libsrc/meshing/python_mesh.cpp

@@ -1480,7 +1480,8 @@ py::arg("point_tolerance") = -1.)
            }))
     */
     
-    .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
+    .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>(),
+          py::arg("dim")=3)
 
     .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
     .def ("BoundaryLayer", [](Mesh & self, variant<string, int, std::vector<int>> boundary,

libsrc/visualization/vssolution.cpp

@@ -841,7 +841,7 @@ namespace netgen
         if (vispar.clipping.enable && clipsolution == 2)      
           {
             mesh->Mutex().unlock();
-            mesh->BuildElementSearchTree();
+            mesh->BuildElementSearchTree(3);
             mesh->Mutex().lock();
           }