From eb8f6f720aaba0446666885fd8154b788e97e32c Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 22 Oct 2015 11:58:13 +0200
Subject: [PATCH] [bboxtree] Use new point-entity intersection detection

Use for pointsources and the collision of a grid with a point
the new algorithms to detect whether a point is inside an entity.
---
 dumux/common/boundingboxtree.hh |  8 ++++++--
 dumux/common/pointsource.hh     | 15 +++++----------
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/dumux/common/boundingboxtree.hh b/dumux/common/boundingboxtree.hh
index 7ce81cce03..19cb58dacc 100644
--- a/dumux/common/boundingboxtree.hh
+++ b/dumux/common/boundingboxtree.hh
@@ -887,11 +887,15 @@ private:
         // We know now it's inside. If the box is a leaf add it.
         else if (isLeaf_(bBox, node))
         {
+            // but add it only if the point is also inside the entity
             const unsigned int eIdx = bBox.child_1;
             auto geometry = (indexToElementMap_->entity(eIdx)).geometry();
-            const ReferenceElement &refElement = ReferenceElements::general(geometry.type());
-            if (refElement.checkInside(geometry.local(point)))
+            if (BoundingBoxTreeHelper<dimworld>::pointInGeometry(geometry, point))
                 entities.push_back(eIdx);
+
+            // const ReferenceElement &refElement = ReferenceElements::general(geometry.type());
+            // if (refElement.checkInside(geometry.local(point)))
+            //     entities.push_back(eIdx);
         }
 
         // No leaf. Check both children.
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 85d1d2b391..e176831ee4 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -25,7 +25,6 @@
 #ifndef DUMUX_POINTSOURCE_HH
 #define DUMUX_POINTSOURCE_HH
 
-#include <dune/geometry/referenceelements.hh>
 #include <dumux/common/boundingboxtree.hh>
 
 namespace Dumux
@@ -79,7 +78,7 @@ private:
 
 /*!
  * \ingroup Common
- * \brief A helper class calculating an DOF index to point source map
+ * \brief A helper class calculating a DOF-index to point source map
  */
 template<class TypeTag>
 class PointSourceHelper
@@ -90,10 +89,7 @@ class PointSourceHelper
     typedef typename GET_PROP_TYPE(TypeTag, PointSource) PointSource;
 
     static const int dim = GridView::dimension;
-
-    typedef typename GridView::Grid::ctype CoordScalar;
-    typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
-    typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
+    static const int dimworld = GridView::dimensionworld;
 
     typedef BoundingBoxTree<GridView> BoundingBoxTree;
 
@@ -119,17 +115,16 @@ public:
                 if(isBox)
                 {
                     // check in which subcontrolvolume(s) we are
-                    auto element = boundingBoxTree->entity(eIdx);
+                    const auto element = boundingBoxTree->entity(eIdx);
                     FVElementGeometry fvGeometry;
                     fvGeometry.update(problem.gridView(), element);
-                    auto globalPos = source.position();
+                    const auto globalPos = source.position();
                     // loop over all sub control volumes and check if the point source is inside
                     std::vector<unsigned int> vertices;
                     for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
                     {
                         auto geometry = fvGeometry.subContVolGeometries[scvIdx];
-                        const ReferenceElement &refElement = ReferenceElements::general(geometry.type());
-                        if (refElement.checkInside(geometry.local(globalPos)))
+                        if (BoundingBoxTreeHelper<dimworld>::pointInGeometry(geometry, globalPos))
                             vertices.push_back(problem.model().dofMapper().subIndex(element, scvIdx, dofCodim));
                     }
                     for (unsigned int vIdx : vertices)
-- 
GitLab