Fix interpolation using rtree.

CMakeLists.txt

@@ -6,7 +6,7 @@ RelWithDebInfo MinSizeRel.")
 PROJECT(Nektar++ C CXX)
 
 CMAKE_POLICY(SET CMP0022 NEW)
-CMAKE_POLICY(SET CMP0074 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")

library/MultiRegions/ExpList.cpp

@@ -1533,28 +1533,38 @@ namespace Nektar
                         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;
-
+//                cout << "Point: " << gloCoords[0] << ", " << gloCoords[1]
+//                     << ", " << gloCoords[2] << endl;
+//                cout << "Number of elements found: " << elmts.size() << endl;
+//                for (int i = 0; i < elmts.size(); ++i) {
+//                    cout << "Element: " << elmts[i] << " of " << GetNumElmts() << endl;
+//                }
                 int min_id  = 0;
                 NekDouble nearpt_min = 1e6;
                 Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
 
                 // search only those elements whose bounding box overlaps
+/*                std::map<int,int> gidToElmtIdxMap;
+                for (int i = 0; i < (*m_exp).size() ; ++i)
+                {
+                    gidToElmtIdxMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
+                }
+*/
                 for (int i = 0; i < elmts.size(); ++i)
                 {
-                    if ((*m_exp)[elmts[i]]->GetGeom()->ContainsPoint(gloCoords, 
+//    cout << "Element Gid: " << elmts[i] << ", index: " << m_elmtToExpId[elmts[i]] << endl;
+                    if ((*m_exp)[m_elmtToExpId[elmts[i]]]->GetGeom()->ContainsPoint(gloCoords, 
                                                               locCoords,
                                                               tol, nearpt))
                     {
-                        return i;
+//                        cout << "LocCoords: " << locCoords[0] << ", " << locCoords[1] << ", " << locCoords[2] << endl;
+                        return m_elmtToExpId[elmts[i]];
                     }
                     else
                     {
                         if(nearpt < nearpt_min)
                         {
-                            min_id    = elmts[i];
+                            min_id    = m_elmtToExpId[elmts[i]];
                             nearpt_min = nearpt;
                             Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
                         }

library/SpatialDomains/Geometry.cpp

@@ -347,20 +347,24 @@ void Geometry::GenBoundingBox()
     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;
+    minx = x - NekConstants::kNekZeroTol;
+    maxx = x + NekConstants::kNekZeroTol;
+    miny = y - NekConstants::kNekZeroTol;
+    maxy = y + NekConstants::kNekZeroTol;
+    minz = z - NekConstants::kNekZeroTol;
+    maxz = z + NekConstants::kNekZeroTol;
     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);
+        miny = (y < miny ? y : miny);
+        maxy = (y > maxy ? y : maxy);
+        minz = (z < minz ? z : minz);
+        maxz = (z > maxz ? z : maxz);
     }
+
     point pmin(minx, miny, minz);
     point pmax(maxx, maxy, maxz);
     m_boundingBox = bg::model::box<point>(pmin, pmax);

library/SpatialDomains/Geometry.h

@@ -49,11 +49,11 @@
 #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
 {
 
+typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> point;
 
 namespace SpatialDomains
 {
@@ -458,7 +458,7 @@ inline bg::model::box<point> Geometry::GetBoundingBox()
 {
     if (!m_boundingBox)
     {
-        GetBoundingBox();
+        GenBoundingBox();
     }
     return *m_boundingBox;
 }

library/SpatialDomains/MeshGraph.cpp

@@ -193,6 +193,11 @@ void MeshGraph::FillBoundingBoxTree()
         box b = x.second->GetBoundingBox();
         m_boundingBoxTree.insert(std::make_pair(b, x.first));
     }
+    for (auto &x : m_quadGeoms) 
+    {
+        box b = x.second->GetBoundingBox();
+        m_boundingBoxTree.insert(std::make_pair(b, x.first));
+    }
 }
 
 void MeshGraph::SetDomainRange(NekDouble xmin, NekDouble xmax, NekDouble ymin,