Initial implementation.

CMakeLists.txt

@@ -6,6 +6,7 @@ RelWithDebInfo MinSizeRel.")
 PROJECT(Nektar++ C CXX)
 
 CMAKE_POLICY(SET CMP0022 NEW)
+CMAKE_POLICY(SET CMP0074 NEW)
 # Nektar++ requires C++11. Try to infer this for older CMake versions (less than
 # 3.1.0)
 IF ("${CMAKE_VERSION}" VERSION_LESS "3.1")

cmake/NektarCommon.cmake

@@ -288,4 +288,4 @@ MACRO(ADD_NEKPY_EXECUTABLE name source)
     INSTALL(FILES ${source} DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
     FILE(COPY ${source} DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
     ADD_CUSTOM_TARGET(${name} SOURCES ${source})
-ENDMACRO()
\ No newline at end of file
+ENDMACRO()

developer-guide @ 0724faa5

-Subproject commit b1461b45abb08c48397fe6d046249703c4f8f160
+Subproject commit 0724faa50ed893acb7ca256b1582d288396bb5d4

library/MultiRegions/ExpList.cpp

@@ -1528,53 +1528,38 @@ namespace Nektar
             // non-embedded mesh (point can only match one element)
             else
             {
-                static int start = 0;
+                SpatialDomains::PointGeomSharedPtr p 
+                    = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(GetExp(0)->GetCoordim(), -1,
+                        gloCoords[0], gloCoords[1], gloCoords[2]);
+
+                std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
+                cout << "Point: " << gloCoords[0] << ", " << gloCoords[1]
+                     << ", " << gloCoords[2] << endl;
+                cout << "Number of elements found: " << elmts.size() << endl;
+
                 int min_id  = 0;
                 NekDouble nearpt_min = 1e6;
                 Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
 
-                // restart search from last found value
-                for (int i = start; i < (*m_exp).size(); ++i)
+                // search only those elements whose bounding box overlaps
+                for (int i = 0; i < elmts.size(); ++i)
                 {
-                    if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, 
+                    if ((*m_exp)[elmts[i]]->GetGeom()->ContainsPoint(gloCoords, 
                                                               locCoords,
                                                               tol, nearpt))
                     {
-                        start = i;
                         return i;
                     }
                     else
                     {
                         if(nearpt < nearpt_min)
                         {
-                            min_id    = i;
+                            min_id    = elmts[i];
                             nearpt_min = nearpt;
                             Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
                         }
                     }
                 }
-
-                for (int i = 0; i < start; ++i)
-                {
-                    if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, 
-                                                              locCoords,
-                                                              tol, nearpt))
-                    {
-                        start = i;
-                        return i;
-                    }
-                    else
-                    {
-                        if(nearpt < nearpt_min)
-                        {
-                            min_id    = i;
-                            nearpt_min = nearpt;
-                            Vmath::Vcopy(locCoords.num_elements(),
-                                         locCoords,1,savLocCoords,1);
-                        }
-                    }
-                }
-
                 if(returnNearestElmt)
                 {
 

library/SpatialDomains/Geometry.cpp

@@ -337,5 +337,36 @@ void Geometry::v_Setup()
              "This function is only valid for expansion type geometries");
 }
 
+
+/**
+ * @brief Generates the bounding box for the element.
+ */
+void Geometry::GenBoundingBox()
+{
+    PointGeomSharedPtr p = GetVertex(0);
+    NekDouble minx, miny, minz, maxx, maxy, maxz;
+    NekDouble x, y, z;
+    p->GetCoords(x, y, z);
+    minx = maxx = x;
+    miny = maxy = y;
+    minz = maxz = z;
+    for (int i = 1; i < GetNumVerts(); ++i)
+    {
+        p = GetVertex(i);
+        p->GetCoords(x, y, z);
+        minx = (x < minx ? x : minx);
+        maxx = (x > maxx ? x : maxx);
+        miny = (y < miny ? x : miny);
+        maxy = (y > maxy ? x : maxy);
+        minz = (z < minz ? x : minz);
+        maxz = (z > maxz ? x : maxz);
+    }
+    point pmin(minx, miny, minz);
+    point pmax(maxx, maxy, maxz);
+    m_boundingBox = bg::model::box<point>(pmin, pmax);
+}
+
+
+
 }
 }

library/SpatialDomains/Geometry.h

@@ -44,9 +44,17 @@
 
 #include <unordered_map>
 
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/geometries/box.hpp>
+namespace bg = boost::geometry;
+namespace bgi = boost::geometry::index;
+typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> point;
+
 namespace Nektar
 {
 
+
 namespace SpatialDomains
 {
 
@@ -138,6 +146,8 @@ public:
     //---------------------------------------
     // Point lookups
     //---------------------------------------
+    SPATIAL_DOMAINS_EXPORT inline bg::model::box<point> GetBoundingBox();
+
     SPATIAL_DOMAINS_EXPORT inline bool ContainsPoint(
         const Array<OneD, const NekDouble> &gloCoord,
         NekDouble tol = 0.0);
@@ -192,8 +202,13 @@ protected:
     int                               m_globalID;
     /// Array containing expansion coefficients of @p m_xmap
     Array<OneD, Array<OneD, NekDouble> > m_coeffs;
+    /// Stores the optional bounding box of the element
+    boost::optional<bg::model::box<point>>   m_boundingBox;
+
     /// Handles generation of geometry factors.
     void                              GenGeomFactors();
+    /// Generates the bounding box of the element.
+    void                              GenBoundingBox();
 
     //---------------------------------------
     // Helper functions
@@ -436,6 +451,18 @@ inline void Geometry::FillGeom()
     v_FillGeom();
 }
 
+/**
+ * @brief Returns the bounding box of the element.
+ */
+inline bg::model::box<point> Geometry::GetBoundingBox()
+{
+    if (!m_boundingBox)
+    {
+        GetBoundingBox();
+    }
+    return *m_boundingBox;
+}
+
 /**
  * @brief Determine whether an element contains a particular Cartesian
  * coordinate \f$(x,y,z)\f$.

library/SpatialDomains/MeshGraph.cpp

@@ -185,6 +185,16 @@ void MeshGraph::FillGraph()
     }
 }
 
+void MeshGraph::FillBoundingBoxTree()
+{
+    m_boundingBoxTree.clear();
+    for (auto &x : m_triGeoms) 
+    {
+        box b = x.second->GetBoundingBox();
+        m_boundingBoxTree.insert(std::make_pair(b, x.first));
+    }
+}
+
 void MeshGraph::SetDomainRange(NekDouble xmin, NekDouble xmax, NekDouble ymin,
                                NekDouble ymax, NekDouble zmin, NekDouble zmax)
 {

library/SpatialDomains/MeshGraph.h

@@ -37,6 +37,10 @@
 #define NEKTAR_SPATIALDOMAINS_MESHGRAPH_H
 
 #include <unordered_map>
+#include <boost/geometry.hpp>
+
+namespace bg = boost::geometry;
+namespace bgi = boost::geometry::index;
 
 #include <LibUtilities/BasicUtils/SessionReader.h>
 #include <LibUtilities/BasicUtils/FieldIO.h>
@@ -157,6 +161,11 @@ typedef std::shared_ptr<std::vector<std::pair<GeometrySharedPtr, int>>>
 
 typedef std::map<std::string, std::string> MeshMetaDataMap;
 
+typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> point;
+typedef bg::model::box<point> box;
+typedef std::pair<box, int> value;
+typedef bgi::rtree< value, bgi::rstar<16, 4> > rtree;
+
 class MeshGraph;
 typedef std::shared_ptr<MeshGraph> MeshGraphSharedPtr;
 
@@ -193,6 +202,28 @@ public:
     /*transfers the minial data structure to full meshgraph*/
     SPATIAL_DOMAINS_EXPORT void FillGraph();
 
+    SPATIAL_DOMAINS_EXPORT void FillBoundingBoxTree();
+    SPATIAL_DOMAINS_EXPORT std::vector<int> GetElementsContainingPoint(
+            PointGeomSharedPtr p)
+    {
+        if (m_boundingBoxTree.empty())
+        {
+            FillBoundingBoxTree();
+        }
+        NekDouble x, y, z;
+        p->GetCoords(x, y, z);
+        box b(point(x, y, z), point(x, y, z));
+        std::vector<value> vals;
+        m_boundingBoxTree.query(bgi::intersects(b), std::back_inserter(vals));
+
+        std::vector<int> ret(vals.size());
+        for (int i = 0; i < ret.size(); ++i)
+        {
+            ret[i] = vals[i].second;
+        }
+        return ret;
+    }
+
     ////////////////////
     ////////////////////
 
@@ -460,6 +491,8 @@ protected:
 
     CompositeOrdering m_compOrder;
     BndRegionOrdering m_bndRegOrder;
+
+    rtree m_boundingBoxTree;
 };
 typedef std::shared_ptr<MeshGraph> MeshGraphSharedPtr;
 typedef LibUtilities::NekFactory<std::string, MeshGraph> MeshGraphFactory;