From e3282300e4ce4ce8bb3249f67f88c59c55a6b872 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 14 Dec 2017 15:49:34 +0100
Subject: [PATCH] [bboxtree] Extract intersection algorithms from bounding box
tree
* Generalize bounding box tree
* Extract intersection algorithms
* Fix some epsilon issues in interval (1d) and triangle (2d) point intersection
* Adapt point sources to new interface
* Adapt fv grid geometry to new interface -> element map moved to basefvgeometry
* Improve bounding box tree test
---
dumux/common/boundingboxtree.hh | 1176 +----------------
dumux/common/elementmap.hh | 50 -
dumux/common/entitymap.hh | 85 ++
dumux/common/geometry/boundingboxtree.hh | 371 ++++++
.../geometry/boundingboxtreeintersection.hh | 93 ++
dumux/common/geometry/geometricentityset.hh | 102 ++
.../geometryintersection.hh} | 88 +-
dumux/common/geometry/intersectingentities.hh | 169 +++
.../geometry/intersectspointgeometry.hh | 111 ++
.../common/geometry/intersectspointsimplex.hh | 249 ++++
dumux/common/pointsource.hh | 16 +-
dumux/discretization/basefvgridgeometry.hh | 34 +-
.../cellcentered/mpfa/fvgridgeometry.hh | 32 +-
.../cellcentered/tpfa/fvgridgeometry.hh | 29 +-
.../staggered/fvgridgeometry.hh | 15 +-
test/common/boundingboxtree/CMakeLists.txt | 14 +-
test/common/boundingboxtree/test_bboxtree.cc | 142 +-
17 files changed, 1375 insertions(+), 1401 deletions(-)
delete mode 100644 dumux/common/elementmap.hh
create mode 100644 dumux/common/entitymap.hh
create mode 100644 dumux/common/geometry/boundingboxtree.hh
create mode 100644 dumux/common/geometry/boundingboxtreeintersection.hh
create mode 100644 dumux/common/geometry/geometricentityset.hh
rename dumux/common/{geometrycollision.hh => geometry/geometryintersection.hh} (65%)
create mode 100644 dumux/common/geometry/intersectingentities.hh
create mode 100644 dumux/common/geometry/intersectspointgeometry.hh
create mode 100644 dumux/common/geometry/intersectspointsimplex.hh
diff --git a/dumux/common/boundingboxtree.hh b/dumux/common/boundingboxtree.hh
index 1b6f24f37d..9c6ed9623b 100644
--- a/dumux/common/boundingboxtree.hh
+++ b/dumux/common/boundingboxtree.hh
@@ -3,7 +3,7 @@
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
+ * the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
@@ -14,1174 +14,8 @@
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
-/*!
- * \file
- * \brief An axis-aligned bounding box volume hierarchy for dune grids
- *
- * Dumux implementation of an AABB tree
- * adapted from implementation in FEniCS by Anders Logg
- */
-#ifndef DUMUX_BOUNDINGBOXTREE_HH
-#define DUMUX_BOUNDINGBOXTREE_HH
-
-#include <algorithm>
-#include <memory>
-#include <type_traits>
-#include <dune/geometry/referenceelements.hh>
-#include <dune/common/timer.hh>
-#include <dune/common/exceptions.hh>
-#include <dumux/common/math.hh>
-#include <dumux/common/exceptions.hh>
-#include <dumux/common/geometrycollision.hh>
-
-namespace Dumux {
-
-/*!
- * \brief A helper class to optimize dimension-dependent methods
- */
-template <int dimworld>
-class BoundingBoxTreeHelper {};
-
-//! Helper methods for world dimension three
-template <>
-class BoundingBoxTreeHelper<3>
-{
- // An epsilon for floating point operations
- static constexpr double eps_ = 1.0e-7;
- using GlobalPosition = Dune::FieldVector<double, 3>;
-
-public:
- //! Check whether a point is inside a given three-dimensional geometry
- template <class Geometry>
- static typename std::enable_if<Geometry::mydimension == 3, bool>::type
- pointInGeometry(const Geometry& geometry, const GlobalPosition& point)
- {
- // get the geometry type
- auto type = geometry.type();
-
- // if it's a tetrahedron we can check directly
- if (type.isTetrahedron())
- {
- return pointInTetrahedron(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(2),
- geometry.corner(3),
- point);
- }
- // split hexahedrons into five tetrahedrons
- else if (type.isHexahedron())
- {
- if (pointInTetrahedron(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(3),
- geometry.corner(5),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(0),
- geometry.corner(5),
- geometry.corner(6),
- geometry.corner(4),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(5),
- geometry.corner(3),
- geometry.corner(6),
- geometry.corner(7),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(0),
- geometry.corner(3),
- geometry.corner(2),
- geometry.corner(6),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(5),
- geometry.corner(3),
- geometry.corner(0),
- geometry.corner(6),
- point))
- return true;
- return false;
- }
- // split pyramids into two tetrahedrons
- else if (type.isPyramid())
- {
- if (pointInTetrahedron(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(2),
- geometry.corner(4),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(1),
- geometry.corner(3),
- geometry.corner(2),
- geometry.corner(4),
- point))
- return true;
- return false;
- }
- // split prisms into three tetrahedrons
- else if (type.isPrism())
- {
- if (pointInTetrahedron(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(2),
- geometry.corner(4),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(3),
- geometry.corner(0),
- geometry.corner(2),
- geometry.corner(4),
- point))
- return true;
- if (pointInTetrahedron(geometry.corner(2),
- geometry.corner(5),
- geometry.corner(3),
- geometry.corner(4),
- point))
- return true;
- return false;
- }
- else
- DUNE_THROW(Dune::NotImplemented,
- "Intersection for point and geometry "
- << type << " in three-dimensional world.");
- }
-
- //! Check whether a point is inside a given two-dimensional geometry
- template <class Geometry>
- static typename std::enable_if<Geometry::mydimension == 2, bool>::type
- pointInGeometry(const Geometry& geometry, const GlobalPosition& point)
- {
- // get the geometry type
- auto type = geometry.type();
-
- // if it's a triangle we can check directly
- if (type.isTriangle())
- {
- return pointInTriangle(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(2),
- point);
- }
- // split quadrilaterals into two triangles
- else if (type.isQuadrilateral())
- {
- if (pointInTriangle(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(3),
- point))
- return true;
- if (pointInTriangle(geometry.corner(0),
- geometry.corner(3),
- geometry.corner(2),
- point))
- return true;
- return false;
- }
- else
- DUNE_THROW(Dune::NotImplemented,
- "Intersection for point and geometry "
- << type << " in three-dimensional world.");
-
- }
-
- //! Check whether a point is inside a given one-dimensional geometry
- template <class Geometry>
- static typename std::enable_if<Geometry::mydimension == 1, bool>::type
- pointInGeometry(const Geometry& geometry, const GlobalPosition& point)
- {
- return pointInInterval(geometry.corner(0),
- geometry.corner(1),
- point);
- }
-
- //! Find out whether a point is inside a tetrahedron (p0, p1, p2, p3)
- static bool pointInTetrahedron(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& p2, const GlobalPosition& p3,
- const GlobalPosition& point)
- {
- // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
- // See also "Real-Time Collision Detection" by Christer Ericson.
- // See also implementation of those algorithms by Anders Logg in FEniCS
-
- // put the tetrahedron points in an array
- const GlobalPosition *p[4] = {&p0, &p1, &p2, &p3};
-
- // iterate over all faces
- for (unsigned int i = 0; i < 4; ++i)
- {
- // compute all the vectors from vertex (local index 0) to the other points
- const GlobalPosition v1 = *p[(i + 1)%4] - *p[i];
- const GlobalPosition v2 = *p[(i + 2)%4] - *p[i];
- const GlobalPosition v3 = *p[(i + 3)%4] - *p[i];
- const GlobalPosition v = point - *p[i];
- // compute the normal to the facet (cross product)
- GlobalPosition n1 = crossProduct(v1, v2);
- n1 /= n1.two_norm();
- // find out on which side of the plane v and v3 are
- const double t1 = n1.dot(v);
- const double t2 = n1.dot(v3);
- // If the point is not exactly on the plane the
- // points have to be on the same side
- const double eps = eps_ * v1.two_norm();
- if ((t1 > eps || t1 < -eps) && std::signbit(t1) != std::signbit(t2))
- return false;
- }
- return true;
- }
-
- //! Find out whether a point is inside a triangle (p0, p1, p2)
- static bool pointInTriangle(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& p2, const GlobalPosition& point)
- {
- // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
- // See also "Real-Time Collision Detection" by Christer Ericson.
- // See also implementation of those algorithms by Anders Logg in FEniCS
-
- // compute the vectors of the edges and the vector point-p0
- const GlobalPosition v1 = p0 - p2;
- const GlobalPosition v2 = p1 - p0;
- const GlobalPosition v3 = p2 - p1;
- const GlobalPosition v = point - p0;
-
- // compute the normal of the triangle
- const GlobalPosition n = crossProduct(v1, v2);
-
- // first check if we are in the plane of the triangle
- // if not we can return early
- const double t = v.dot(n);
- using std::abs;
- if (abs(t) > v1.two_norm()*eps_) // take |v1| as scale
- return false;
-
- // compute the normal to the triangle made of point and first edge
- // the dot product of this normal and the triangle normal has to
- // be positive because we defined the edges in the right orientation
- const GlobalPosition n1 = crossProduct(v, v1);
- const double t1 = n.dot(n1);
- if (t1 < 0) return false;
-
- const GlobalPosition n2 = crossProduct(v, v2);
- const double t2 = n.dot(n2);
- if (t2 < 0) return false;
-
- const GlobalPosition n3 = crossProduct(v, v3);
- const double t3 = n.dot(n3);
- if (t3 < 0) return false;
-
- return true;
- }
-
- //! Find out whether a point is inside an interval (p0, p1)
- static bool pointInInterval(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& point)
- {
- // compute the vectors between p0 and the other points
- const GlobalPosition v1 = p1 - p0;
- const GlobalPosition v2 = point - p0;
-
- // check if point and p0 are the same
- const double v1norm = v1.two_norm();
- const double v2norm = v2.two_norm();
- if (v2norm < v1norm*eps_)
- return true;
-
- // if not check if p0 and p1 are the same
- // then we know that point is not in the interval
- if (v1norm < eps_)
- return false;
-
- // if the cross product is zero the points are on a line
- const GlobalPosition n = crossProduct(v1, v2);
-
- // early return if the vector length is larger than zero
- if (n.two_norm() > v1norm*eps_)
- return false;
-
- // we know the points are aligned
- // if the dot product is positive and the length in range
- // the point is in the interval
- return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
- }
-
- /*!
- * \brief Check whether a point is in a bounding box
- * \param point The point
- * \param b Pointer to bounding box coordinates
- */
- static bool pointInBoundingBox(const Dune::FieldVector<double, 3>& point,
- const double* b)
- {
- const double eps0 = eps_*(b[3] - b[0]);
- const double eps1 = eps_*(b[4] - b[1]);
- const double eps2 = eps_*(b[5] - b[2]);
- return (b[0] - eps0 <= point[0] && point[0] <= b[3] + eps0 &&
- b[1] - eps1 <= point[1] && point[1] <= b[4] + eps1 &&
- b[2] - eps2 <= point[2] && point[2] <= b[5] + eps2);
- }
-
- /*!
- * \brief Check whether bounding box a collides with bounding box b
- * \param a, b Pointer to bounding box coordinates
- */
- static bool boundingBoxInBoundingBox(const double* a,
- const double* b)
- {
- const double eps0 = eps_*(b[3] - b[0]);
- const double eps1 = eps_*(b[4] - b[1]);
- const double eps2 = eps_*(b[5] - b[2]);
- return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
- b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
- b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
- }
-
- //! Compute the bounding box of a vector of bounding boxes
- static void computeBBoxOfBBoxes(double* bBox,
- std::size_t& axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& end)
- {
- // Copy the iterator and get coordinates of first box
- auto it = begin;
- const double* bFirst = leafBoxes.data() + 6*(*it);
- // Maybe write out the loop for optimization
- for (std::size_t coordIdx = 0; coordIdx < 6; ++coordIdx)
- bBox[coordIdx] = bFirst[coordIdx];
-
- // Compute min and max with the remaining boxes
- for (; it != end; ++it)
- {
- const double* b = leafBoxes.data() + 6*(*it);
- if (b[0] < bBox[0]) bBox[0] = b[0];
- if (b[1] < bBox[1]) bBox[1] = b[1];
- if (b[2] < bBox[2]) bBox[2] = b[2];
- if (b[3] > bBox[3]) bBox[3] = b[3];
- if (b[4] > bBox[4]) bBox[4] = b[4];
- if (b[5] > bBox[5]) bBox[5] = b[5];
- }
-
- // Compute the longest axis
- const double x = bBox[3] - bBox[0];
- const double y = bBox[4] - bBox[1];
- const double z = bBox[5] - bBox[2];
-
- if (x > y && x > z)
- axis = 0;
- else if (y > z)
- axis = 1;
- else
- axis = 2;
- }
-
- //! Sort the bounding boxes along the longest axis
- static void sortBoundingBoxes(std::size_t axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& middle,
- const std::vector<unsigned int>::iterator& end)
- {
- switch (axis)
- {
- case 0:
- std::nth_element(begin, middle, end, lessXBox(leafBoxes));
- case 1:
- std::nth_element(begin, middle, end, lessYBox(leafBoxes));
- default:
- std::nth_element(begin, middle, end, lessZBox(leafBoxes));
- }
- }
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the x-axis
- * \note This could be replaced by lambdas
- */
- struct lessXBox
- {
- const std::vector<double>& bBoxes;
- lessXBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 6*i;
- const double* bj = bBoxes.data() + 6*j;
- return bi[0] + bi[3] < bj[0] + bj[3];
- }
- };
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the y-axis
- * \note This could be replaced by lambdas
- */
- struct lessYBox
- {
- const std::vector<double>& bBoxes;
- lessYBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 6*i;
- const double* bj = bBoxes.data() + 6*j;
- return bi[1] + bi[4] < bj[1] + bj[4];
- }
- };
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the z-axis
- * \note This could be replaced by lambdas
- */
- struct lessZBox
- {
- const std::vector<double>& bBoxes;
- lessZBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 6*i;
- const double* bj = bBoxes.data() + 6*j;
- return bi[2] + bi[5] < bj[2] + bj[5];
- }
- };
-};
-
-//! Helper methods for world dimension two
-template <>
-class BoundingBoxTreeHelper<2>
-{
- // An epsilon for floating point operations
- static constexpr double eps_ = 1.0e-7;
- using GlobalPosition = Dune::FieldVector<double, 2>;
-
-public:
- // Check whether a point is inside a given geometry
- template <class Geometry>
- static typename std::enable_if<Geometry::mydimension == 2, bool>::type
- pointInGeometry(const Geometry& geometry, const GlobalPosition& point)
- {
- // get the geometry type
- auto type = geometry.type();
-
- // if it's a triangle we can check directly
- if (type.isTriangle())
- {
- return pointInTriangle(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(2),
- point);
- }
- // split quadrilaterals into two triangles
- else if (type.isQuadrilateral())
- {
- if (pointInTriangle(geometry.corner(0),
- geometry.corner(1),
- geometry.corner(3),
- point))
- return true;
- if (pointInTriangle(geometry.corner(0),
- geometry.corner(3),
- geometry.corner(2),
- point))
- return true;
- return false;
- }
- else
- DUNE_THROW(Dune::NotImplemented,
- "Intersection for point and geometry "
- << type << " in two-dimensional world.");
- }
-
- template <class Geometry>
- static typename std::enable_if<Geometry::mydimension == 1, bool>::type
- pointInGeometry(const Geometry& geometry, const GlobalPosition& point)
- {
- return pointInInterval(geometry.corner(0),
- geometry.corner(1),
- point);
- }
-
- //! find out if a point is inside a triangle
- static bool pointInTriangle(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& p2, const GlobalPosition& point)
- {
- // Use barycentric coordinates
- const double A = 0.5*(-p1[1]*p2[0] + p0[1]*(p2[0] - p1[0])
- + p0[0]*(p1[1] - p2[1]) + p1[0]*p2[1]);
- const double sign = std::copysign(1.0, A);
- const double s = sign*(p0[1]*p2[0] + point[0]*(p2[1]-p0[1])
- -p0[0]*p2[1] + point[1]*(p0[0]-p2[0]));
- const double t = sign*(p0[0]*p1[1] + point[0]*(p0[1]-p1[1])
- -p0[1]*p1[0] + point[1]*(p1[0]-p0[0]));
- const double eps = (p0 - p1).two_norm()*eps_;
-
- return (s > -eps
- && t > -eps
- && (s + t) < 2*A*sign + eps);
- }
-
- //! find out if a point is inside an interval
- static bool pointInInterval(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& point)
- {
- // compute the vectors between p0 and the other points
- const GlobalPosition v1 = p1 - p0;
- const GlobalPosition v2 = point - p0;
-
- // check if point and p0 are the same
- const double v1norm = v1.two_norm();
- const double v2norm = v2.two_norm();
- if (v2norm < v1norm*eps_)
- return true;
-
- // if not check if p0 and p1 are the same
- // then we know that point is not in the interval
- if (v1norm < eps_)
- return false;
-
- // if the cross product is zero the points are on a line
- const double n = crossProduct(v1, v2);
-
- // early return if the cross product is larger than zero
- using std::abs;
- if (abs(n) > v1norm*eps_)
- return false;
-
- // we know the points are aligned
- // if the dot product is positive and the length in range
- // the point is in the interval
- return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
- }
-
- /*!
- * \brief Check whether a point is in a bounding box
- * \param point The point
- * \param b Pointer to bounding box coordinates
- */
- static bool pointInBoundingBox(const Dune::FieldVector<double, 2>& point,
- const double* b)
- {
- const double eps0 = eps_*(b[2] - b[0]);
- const double eps1 = eps_*(b[3] - b[1]);
- return (b[0] - eps0 <= point[0] && point[0] <= b[2] + eps0 &&
- b[1] - eps1 <= point[1] && point[1] <= b[3] + eps1);
- }
-
- /*!
- * \brief Check whether bounding box a collides with bounding box b
- * \param a, b Pointer to bounding box coordinates
- */
- static bool boundingBoxInBoundingBox(const double* a,
- const double* b)
- {
- const double eps0 = eps_*(b[2] - b[0]);
- const double eps1 = eps_*(b[3] - b[1]);
- return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
- b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
- }
-
- //! Compute the bounding box of a vector of bounding boxes
- static void computeBBoxOfBBoxes(double* bBox,
- std::size_t& axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& end)
- {
- // Copy the iterator and get coordinates of first box
- auto it = begin;
- const double* bFirst = leafBoxes.data() + 4*(*it);
- // Maybe write out the loop for optimization
- for (std::size_t coordIdx = 0; coordIdx < 4; ++coordIdx)
- bBox[coordIdx] = bFirst[coordIdx];
-
- // Compute min and max with the remaining boxes
- for (; it != end; ++it)
- {
- const double* b = leafBoxes.data() + 4*(*it);
- if (b[0] < bBox[0]) bBox[0] = b[0];
- if (b[1] < bBox[1]) bBox[1] = b[1];
- if (b[2] > bBox[2]) bBox[2] = b[2];
- if (b[3] > bBox[3]) bBox[3] = b[3];
- }
-
- // Compute the longest axis
- const double x = bBox[2] - bBox[0];
- const double y = bBox[3] - bBox[1];
-
- if (x > y)
- axis = 0;
- else
- axis = 1;
- }
-
- //! Sort the bounding boxes along the longest axis
- static void sortBoundingBoxes(std::size_t axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& middle,
- const std::vector<unsigned int>::iterator& end)
- {
- if (axis == 0)
- std::nth_element(begin, middle, end, lessXBox(leafBoxes));
- else
- std::nth_element(begin, middle, end, lessYBox(leafBoxes));
- }
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the x-axis
- * \note This could be replaced by lambdas
- */
- struct lessXBox
- {
- const std::vector<double>& bBoxes;
- lessXBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 4*i;
- const double* bj = bBoxes.data() + 4*j;
- return bi[0] + bi[2] < bj[0] + bj[2];
- }
- };
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the y-axis
- * \note This could be replaced by lambdas
- */
- struct lessYBox
- {
- const std::vector<double>& bBoxes;
- lessYBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 4*i;
- const double* bj = bBoxes.data() + 4*j;
- return bi[1] + bi[3] < bj[1] + bj[3];
- }
- };
-};
-
-//! Helper methods for world dimension one
-template <>
-class BoundingBoxTreeHelper<1>
-{
- // An epsilon for floating point operations
- static constexpr double eps_ = 1.0e-7;
- using GlobalPosition = Dune::FieldVector<double, 1>;
-
-public:
- // Check whether a point is inside a given geometry
- template <class Geometry>
- static bool pointInGeometry(const Geometry& geometry,
- const GlobalPosition& point)
- {
- return pointInInterval(geometry.corner(0),
- geometry.corner(1),
- point);
- }
-
- //! find out if a point is inside an interval
- static bool pointInInterval(const GlobalPosition& p0, const GlobalPosition& p1,
- const GlobalPosition& point)
- {
- const double v1 = point[0] - p0[0];
- const double v2 = p1[0] - p0[0];
- // the coordinates are the same
- if (v1 < v2*eps_)
- return true;
- // the point doesn't coincide with p0
- // so if p0 and p1 are equal it's not inside
- if (v2 < eps_)
- return false;
-
- // the point is inside if the length is
- // small than the interval length and the
- // sign of v1 & v2 are the same
- using std::abs;
- using std::signbit;
- return (signbit(v1) == signbit(v2)
- && abs(v1) < abs(v2)*(1 + eps_));
- }
-
- /*!
- * \brief Check whether a point is in a bounding box
- * \param point The point
- * \param b Pointer to bounding box coordinates
- */
- static bool pointInBoundingBox(const GlobalPosition& point,
- const double* b)
- {
- const double eps0 = eps_*(b[1] - b[0]);
- return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
- }
-
- /*!
- * \brief Check whether bounding box a collides with bounding box b
- * \param a, b Pointer to bounding box coordinates
- */
- static bool boundingBoxInBoundingBox(const double* a,
- const double* b)
- {
- const double eps0 = eps_*(b[1] - b[0]);
- return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
- }
-
- //! Compute the bounding box of a vector of bounding boxes
- static void computeBBoxOfBBoxes(double* bBox,
- std::size_t& axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& end)
- {
- // Copy the iterator and get coordinates of first box
- auto it = begin;
- const double* bFirst = leafBoxes.data() + 2*(*it);
- bBox[0] = bFirst[0];
- bBox[1] = bFirst[1];
-
- // Compute min and max with the remaining boxes
- for (; it != end; ++it)
- {
- const double* b = leafBoxes.data() + 2*(*it);
- if (b[0] < bBox[0]) bBox[0] = b[0];
- if (b[1] > bBox[1]) bBox[1] = b[1];
- }
-
- // Compute the longest axis
- axis = 0;
- }
-
- //! Sort the bounding boxes along the longest axis
- static void sortBoundingBoxes(std::size_t axis,
- const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& middle,
- const std::vector<unsigned int>::iterator& end)
- { std::nth_element(begin, middle, end, lessXBox(leafBoxes)); }
-
- /*!
- * \brief Comparison function for sorting bounding boxes on the x-axis
- * \note This could be replaced by lambdas
- */
- struct lessXBox
- {
- const std::vector<double>& bBoxes;
- lessXBox(const std::vector<double>& bBoxes_): bBoxes(bBoxes_) {}
- inline bool operator()(unsigned int i, unsigned int j)
- {
- const double* bi = bBoxes.data() + 2*i;
- const double* bj = bBoxes.data() + 2*j;
- return bi[0] + bi[1] < bj[0] + bj[1];
- }
- };
-};
-
-/*!
- * \brief An intersection object resulting from the collision of two bounding box tree primitives
- *
- * After is has been found that two leaf bounding boxes intersect a primitive test has to be
- * performed to see if the actual entities inside the bounding box intersect too. The result
- * if such an intersection is found as an object of this class containing the indices of the
- * intersecting entities and the corners of the intersection object.
- */
-template<class GridView1, class GridView2>
-class BoundingBoxTreeIntersection
-{
- const static int dimworld = GridView1::dimensionworld;
- using Scalar = typename GridView1::ctype;
- using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
-
-public:
- BoundingBoxTreeIntersection(unsigned int a, unsigned int b, std::vector<GlobalPosition>&& c)
- : a_(a), b_(b), corners_(c) {}
-
- //! Get the index of the intersecting entity belonging to this grid
- inline unsigned int first() const
- { return a_; }
-
- //! Get the index of the intersecting entity belonging to the other grid
- inline unsigned int second() const
- { return b_; }
-
- //! Get the corners of the intersection geometry
- inline std::vector<GlobalPosition> corners() const
- { return corners_; }
-
- /*!
- * \brief Check if the corners of this intersection match with the given corners
- * \note This is useful to check if the intersection geometry of two intersections coincide.
- */
- bool cornersMatch(const std::vector<GlobalPosition>& corners) const
- {
- if (corners.size() != corners_.size())
- return false;
-
- const auto eps = 1.5e-7*(corners_[1] - corners_[0]).two_norm();
- for (int i = 0; i < corners_.size(); ++i)
- if ((corners[0] - corners[1]).two_norm() > eps)
- return false;
-
- return true;
- }
-
-private:
- unsigned int a_, b_; //!< Indices of the intersection elements
- std::vector<GlobalPosition> corners_; //!< the corner points of the intersection geometry
-};
-
-/*!
- * \brief A class mapping from an element index to elements using element seeds
- */
-template <class GridView>
-class IndexToElementMap
- : public std::vector<typename GridView::Traits::Grid::template Codim<0>::EntitySeed>
-{
- using Grid = typename GridView::Traits::Grid;
- using Element = typename GridView::template Codim<0>::Entity;
-public:
- IndexToElementMap(const GridView& gridView)
- : grid_(gridView.grid()) {}
-
- //! get an element from an index t
- template<class T>
- Element entity(T&& t)
- { return grid_.entity((*this)[std::forward<T>(t)]); }
-
-private:
- const Grid& grid_;
-};
-
-/*!
- * \brief An axis-aligned bounding box volume tree implementation
- *
- * The class constructs a hierarchical structure of bounding box volumes around
- * grid entities. This class can be used to efficiently compute intersections
- * between a grid and other geometrical object. It only implements the intersection
- * of two of such bounding box trees, so that two independent grids can be intersected.
- */
-template <class GridView>
-class BoundingBoxTree
-{
- //! be friends with all other kinds bounding box trees so that
- //! they can call each others private methods
- template <class OtherGridView> friend class BoundingBoxTree;
-
- static const int dim = GridView::dimension;
- static const int dimworld = GridView::dimensionworld;
- using Element = typename GridView::template Codim<0>::Entity;
- using Scalar = typename GridView::ctype;
- using GlobalPosition = typename Dune::FieldVector<Scalar, dimworld>;
-
-public:
- //! Default Constructor
- BoundingBoxTree() {}
-
- //! Constructor with gridView
- BoundingBoxTree(const GridView& leafGridView)
- { build(leafGridView); }
-
- //! Build up bounding box tree for a grid with leafGridView
- void build(const GridView& leafGridView)
- {
- // clear data if any
- clear_();
-
- // Create bounding boxes for all elements
- const unsigned int numLeaves = leafGridView.size(0);
-
- // start the timer
- Dune::Timer timer;
-
- // create a vector for leaf boxes (min and max for all dims)
- std::vector<double> leafBoxes(2*dimworld*numLeaves);
-
- // Create index to element map
- indexToElementMap_ = std::make_shared<IndexToElementMap<GridView> >(leafGridView);
- indexToElementMap_->resize(numLeaves);
-
- for (auto&& element : elements(leafGridView))
- {
- unsigned int eIdx = leafGridView.indexSet().index(element);
- computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*eIdx, element);
- (*indexToElementMap_)[eIdx] = element.seed();
- }
-
- // Create the leaf partition, the set of available indices (to be sorted)
- std::vector<unsigned int> leafPartition(numLeaves);
- for (unsigned int i = 0; i < numLeaves; ++i)
- leafPartition[i] = i;
-
- // Recursively build the bounding box tree
- build_(leafBoxes, leafPartition.begin(), leafPartition.end());
-
- // We are done, log output
- std::cout << "Computed bounding box tree with " << numBoundingBoxes_()
- << " nodes for " << numLeaves << " grid entites in "
- << timer.stop() << " seconds." << std::endl;
- }
-
- //! Compute all intersections between entities and a point
- std::vector<unsigned int> computeEntityCollisions(const Dune::FieldVector<double, dimworld>& point) const
- {
- // Call the recursive find function to find candidates
- std::vector<unsigned int> entities;
- computeCollisions_(point, numBoundingBoxes_() - 1, entities);
- return entities;
- }
-
- //! Compute all intersections between entities and another bounding box tree
- template<class OtherGridView>
- std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>>
- computeEntityCollisions(const BoundingBoxTree<OtherGridView>& otherTree) const
- {
- // check if the world dimensions match
- static_assert(dimworld == OtherGridView::dimensionworld, "Can only collide bounding box trees of same world dimension");
-
- // Create data structure for return type
- std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>> intersections;
-
- // Call the recursive find function to find candidates
- computeCollisions_(otherTree,
- this->numBoundingBoxes_() - 1,
- otherTree.numBoundingBoxes_() -1,
- intersections);
-
- return intersections;
- }
-
- //! Get an element from a global element index
- Element entity(unsigned int eIdx) const
- { return indexToElementMap_->entity(eIdx); }
-
-private:
- /*!
- * \brief Bounding box data structure
- * Leaf nodes are indicated by setting child_0 to
- * the node itself and child_1 is the index of the entity in the bounding box.
- */
- struct BoundingBox
- {
- unsigned int child_0;
- unsigned int child_1;
- };
-
- //! Vector of bounding boxes
- std::vector<BoundingBox> boundingBoxes_;
-
- //! Vector of bounding box coordinates
- std::vector<double> boundingBoxCoordinates_;
-
- //! Shared pointer to the index to element map
- std::shared_ptr<IndexToElementMap<GridView> > indexToElementMap_;
-
- //! Clear all data
- void clear_()
- {
- boundingBoxes_.clear();
- boundingBoxCoordinates_.clear();
- if(indexToElementMap_) indexToElementMap_->clear();
- }
-
- //! Build bounding box tree for all entities recursively
- unsigned int build_(const std::vector<double>& leafBoxes,
- const std::vector<unsigned int>::iterator& begin,
- const std::vector<unsigned int>::iterator& end)
- {
- assert(begin < end);
-
- // Create empty bounding box
- BoundingBox bBox;
-
- // If we reached the leaf
- if (end - begin == 1)
- {
- // Get the bounding box coordinates for the leaf
- const unsigned int eIdx = *begin;
- const double* b = leafBoxes.data() + 2*dimworld*eIdx;
-
- // Store the data in the bounding box
- bBox.child_0 = numBoundingBoxes_();
- bBox.child_1 = eIdx;
- return addBoundingBox_(bBox, b);
- }
-
- // Compute the bounding box of all bounding boxes in the range [begin, end]
- double b[dimworld*2];
- std::size_t axis;
- BoundingBoxTreeHelper<dimworld>::computeBBoxOfBBoxes(b, axis, leafBoxes, begin, end);
-
- // Sort bounding boxes along the longest axis
- std::vector<unsigned int>::iterator middle = begin + (end - begin)/2;
- BoundingBoxTreeHelper<dimworld>::sortBoundingBoxes(axis, leafBoxes, begin, middle, end);
-
- // Split the bounding boxes into two branches and call build recursively
- bBox.child_0 = build_(leafBoxes, begin, middle);
- bBox.child_1 = build_(leafBoxes, middle, end);
-
- // Store the bounding box data. The root will be added at the end.
- return addBoundingBox_(bBox, b);
- }
-
- //! Compute collisions with point recursively
- void computeCollisions_(const Dune::FieldVector<double, dimworld>& point,
- unsigned int node,
- std::vector<unsigned int>& entities) const
- {
- // Get the bounding box for the current node
- const BoundingBox& bBox = getBoundingBox_(node);
-
- // if the point is not in the bounding box we can stop
- if (!BoundingBoxTreeHelper<dimworld>::pointInBoundingBox(point, getBoundingBoxCoordinates_(node)))
- return;
-
- // 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();
- 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.
- else
- {
- computeCollisions_(point, bBox.child_0, entities);
- computeCollisions_(point, bBox.child_1, entities);
- }
- }
-
- //! Compute collisions with other bounding box tree recursively
- template <class OtherGridView>
- void computeCollisions_(const BoundingBoxTree<OtherGridView>& treeB,
- unsigned int nodeA,
- unsigned int nodeB,
- std::vector<BoundingBoxTreeIntersection<GridView, OtherGridView>>& intersections) const
- {
- // get alias
- const auto& treeA = *this;
-
- // Get the bounding box for the current node
- const auto& bBoxA = treeA.getBoundingBox_(nodeA);
- const auto& bBoxB = treeB.getBoundingBox_(nodeB);
-
- // if the two bounding boxes don't collide we can stop searching
- if (!BoundingBoxTreeHelper<dimworld>::
- boundingBoxInBoundingBox(treeA.getBoundingBoxCoordinates_(nodeA),
- treeB.getBoundingBoxCoordinates_(nodeB)))
- return;
-
- // Check if we have a leaf in treeA or treeB
- const bool isLeafA = treeA.isLeaf_(bBoxA, nodeA);
- const bool isLeafB = treeB.isLeaf_(bBoxB, nodeB);
-
- // If both boxes are leaves and collide add them
- if (isLeafA && isLeafB)
- {
- const unsigned int eIdxA = bBoxA.child_1;
- const unsigned int eIdxB = bBoxB.child_1;
-
- auto geometryA = treeA.entity(eIdxA).geometry();
- auto geometryB = treeB.entity(eIdxB).geometry();
-
- using CollisionType = GeometryCollision<decltype(geometryA), decltype(geometryB)>;
- std::vector<GlobalPosition> intersection;
- if (CollisionType::collide(geometryA, geometryB, intersection))
- intersections.emplace_back(eIdxA, eIdxB, std::move(intersection));
- }
-
- // if we reached the leaf in treeA, just continue in treeB
- else if (isLeafA)
- {
- computeCollisions_(treeB, nodeA, bBoxB.child_0, intersections);
- computeCollisions_(treeB, nodeA, bBoxB.child_1, intersections);
- }
-
- // if we reached the leaf in treeB, just continue in treeA
- else if (isLeafB)
- {
- computeCollisions_(treeB, bBoxA.child_0, nodeB, intersections);
- computeCollisions_(treeB, bBoxA.child_1, nodeB, intersections);
- }
-
- // we know now that both trees didn't reach the leaf yet so
- // we continue with the larger tree first (bigger node number)
- else if (nodeA > nodeB)
- {
- computeCollisions_(treeB, bBoxA.child_0, nodeB, intersections);
- computeCollisions_(treeB, bBoxA.child_1, nodeB, intersections);
- }
- else
- {
- computeCollisions_(treeB, nodeA, bBoxB.child_0, intersections);
- computeCollisions_(treeB, nodeA, bBoxB.child_1, intersections);
- }
- }
-
- //! Add a new bounding box to the tree
- inline unsigned int addBoundingBox_(const BoundingBox& bBox,
- const double* b)
- {
- // Add the bounding box
- boundingBoxes_.push_back(bBox);
-
- // Add the bounding box's coordinates
- for (std::size_t i = 0; i < 2*dimworld; ++i)
- boundingBoxCoordinates_.push_back(b[i]);
-
- // return the index of the new node
- return boundingBoxes_.size() - 1;
- }
-
- //! Get an existing bounding box for a given node
- inline const BoundingBox& getBoundingBox_(unsigned int node) const
- { return boundingBoxes_[node]; }
-
- //! Get an existing bounding box for a given node
- const double* getBoundingBoxCoordinates_(unsigned int node) const
- { return boundingBoxCoordinates_.data() + 2*dimworld*node; }
-
- //! Get the number of bounding boxes currently in the tree
- inline std::size_t numBoundingBoxes_() const
- { return boundingBoxes_.size(); }
-
- //! Check whether a bounding box is a leaf node
- //! Leaf nodes have itself as child_0
- inline bool isLeaf_(const BoundingBox& bBox, unsigned int node) const
- { return bBox.child_0 == node; }
-
- //! Compute the bounding box of a grid entity
- template <class Entity>
- void computeEntityBoundingBox_(double* b, const Entity& entity) const
- {
- // get the bounding box coordinates
- double* xMin = b;
- double* xMax = b + dimworld;
-
- // get mesh entity data
- auto geometry = entity.geometry();
-
- // Get coordinates of first vertex
- auto corner = geometry.corner(0);
- for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
- xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
-
- // Compute the min and max over the remaining vertices
- for (std::size_t vLocalIdx = 1; vLocalIdx < entity.subEntities(dim); ++vLocalIdx)
- {
- corner = geometry.corner(vLocalIdx);
- for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
- {
- using std::max;
- using std::min;
- xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
- xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
- }
- }
- }
-};
-
-} // end namespace Dumux
-
+#ifndef DUMUX_OLD_BOUNDINGBOXTREE_HH
+#define DUMUX_OLD_BOUNDINGBOXTREE_HH
+#warning "This header is deprecated. Include dumux/common/geometry/boundingboxtree.hh"
+#include <dumux/common/geometry/boundingboxtree.hh>
#endif
diff --git a/dumux/common/elementmap.hh b/dumux/common/elementmap.hh
deleted file mode 100644
index 5ffb0c26b6..0000000000
--- a/dumux/common/elementmap.hh
+++ /dev/null
@@ -1,50 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Base class for the finite volume geometry vector for cell-centered TPFA models
- * This builds up the sub control volumes and sub control volume faces
- * for each element.
- */
-#ifndef DUMUX_ELEMENT_MAP_HH
-#define DUMUX_ELEMENT_MAP_HH
-
-namespace Dumux {
-
-//! An index to element map
-template <class GridView>
-class ElementMap
- : public std::vector<typename GridView::Traits::Grid::template Codim<0>::EntitySeed>
-{
- using Grid = typename GridView::Traits::Grid;
- using Element = typename GridView::template Codim<0>::Entity;
- using IndexType = typename GridView::IndexSet::IndexType;
-public:
- ElementMap(const GridView& gridView_) : grid_(gridView_.grid()) {}
-
- Element element(IndexType eIdx) const
- { return grid_.entity((*this)[eIdx]); }
-
-private:
- const Grid& grid_;
-};
-
-} // end namespace Dumux
-
-#endif
\ No newline at end of file
diff --git a/dumux/common/entitymap.hh b/dumux/common/entitymap.hh
new file mode 100644
index 0000000000..de4a42b773
--- /dev/null
+++ b/dumux/common/entitymap.hh
@@ -0,0 +1,85 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A map from indices to entities using grid entity seeds
+ */
+#ifndef DUMUX_ENTITY_INDEX_MAP_HH
+#define DUMUX_ENTITY_INDEX_MAP_HH
+
+namespace Dumux {
+
+/*!
+ * \brief A map from indices to entities using grid entity seeds
+ */
+template <class GridView, int codim = 0>
+class EntityMap
+{
+public:
+ using Grid = typename GridView::Traits::Grid;
+ using Entity = typename Grid::template Codim<codim>::Entity;
+ using EntitySeed = typename Grid::template Codim<codim>::EntitySeed;
+
+ //! constructor moving a ready seed list in here
+ EntityMap(const Grid& grid, std::vector<EntitySeed>&& seeds)
+ : grid_(grid)
+ , seeds_(std::move(seeds))
+ {}
+
+ //! constructor with all entites of codim
+ template<class Mapper>
+ EntityMap(const Grid& grid, const Mapper& mapper)
+ : grid_(grid)
+ {
+ update(mapper);
+ }
+
+ //! update the map after the grid changed
+ void update(std::vector<EntitySeed>&& seeds)
+ { seeds_.swap(std::move(seeds)); }
+
+ //! update the map after the grid changed
+ template<class Mapper>
+ void update(const Mapper& mapper)
+ {
+ const auto& gv = grid_.leafGridView();
+ seeds_.resize(gv.size(codim));
+ for (const auto& entity : entities(gv, Dune::Codim<codim>()))
+ seeds_[mapper.index(entity)] = entity.seed();
+ }
+
+ //! get an element from an index i
+ Entity operator[](std::size_t i) const
+ { return grid_.entity(seeds_[i]); }
+
+ //! get the size of the map
+ std::size_t size() const
+ { return seeds_.size(); }
+
+private:
+ const Grid& grid_;
+ std::vector<EntitySeed> seeds_;
+};
+
+template<class GridView>
+using ElementMap = EntityMap<GridView, 0>;
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/boundingboxtree.hh b/dumux/common/geometry/boundingboxtree.hh
new file mode 100644
index 0000000000..e3aec399d1
--- /dev/null
+++ b/dumux/common/geometry/boundingboxtree.hh
@@ -0,0 +1,371 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An axis-aligned bounding box volume hierarchy for dune grids
+ *
+ * Dumux implementation of an AABB tree
+ * Inspired by the AABB tree implementation in DOLFIN by Anders Logg which has the
+ * following license info: DOLFIN is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ */
+#ifndef DUMUX_BOUNDINGBOXTREE_HH
+#define DUMUX_BOUNDINGBOXTREE_HH
+
+#include <vector>
+#include <array>
+#include <algorithm>
+#include <numeric>
+#include <type_traits>
+
+#include <dune/common/promotiontraits.hh>
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \brief An axis-aligned bounding box volume tree implementation
+ *
+ * The class constructs a hierarchical structure of bounding box volumes around
+ * grid entities. This class can be used to efficiently compute intersections
+ * between a grid and other geometrical object. It only implements the intersection
+ * of two of such bounding box trees, so that two independent grids can be intersected.
+ * \tparam GeometricEntitySet has the following requirements
+ * * export dimensionworld, ctype
+ * * a size() member function returning the number of entities
+ * * begin() and end() member function returning at least forward iterators to entities
+ * * an index() method returning a consecutive index given an entity
+ * * an entity() method returning an entity given the consecutive index
+ * * entities have the following requirements:
+ * * a member function geometry() returning a geometry with the member functions
+ * * corner() and corners() returning global coordinates and number of corners
+ */
+template <class GeometricEntitySet>
+class BoundingBoxTree
+{
+ enum { dimworld = GeometricEntitySet::dimensionworld };
+ using ctype = typename GeometricEntitySet::ctype;
+
+ /*!
+ * \brief Bounding box node data structure
+ * leaf nodes are indicated by setting child0 to
+ * the node itself and child1 to the index of the entity in the bounding box.
+ */
+ struct BoundingBoxNode
+ {
+ std::size_t child0;
+ std::size_t child1;
+ };
+
+public:
+ //! the type of entity set this tree was built with
+ using EntitySet = GeometricEntitySet;
+
+ //! Default Constructor
+ BoundingBoxTree() = default;
+
+ //! Constructor with gridView
+ BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set)
+ { build(set); }
+
+ //! Build up bounding box tree for a grid with leafGridView
+ void build(std::shared_ptr<const GeometricEntitySet> set)
+ {
+ // set the pointer to the entity set
+ entitySet_ = set;
+
+ // clear all internal data
+ boundingBoxNodes_.clear();
+ boundingBoxCoordinates_.clear();
+
+ // start the timer
+ Dune::Timer timer;
+
+ // Create bounding boxes for all elements
+ const auto numLeaves = set->size();
+
+ // reserve enough space for the nodes and the coordinates
+ const auto numNodes = 2*numLeaves - 1;
+ boundingBoxNodes_.reserve(numNodes);
+ boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
+
+ // create a vector for leaf boxes (min and max for all dims)
+ std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
+
+ for (const auto& geometricEntity : *set)
+ computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
+
+ // create the leaf partition, the set of available indices (to be sorted)
+ std::vector<std::size_t> leafPartition(numLeaves);
+ std::iota(leafPartition.begin(), leafPartition.end(), 0);
+
+ // Recursively build the bounding box tree
+ build_(leafBoxes, leafPartition.begin(), leafPartition.end());
+
+ // We are done, log output
+ std::cout << "Computed bounding box tree with " << numBoundingBoxes()
+ << " nodes for " << numLeaves << " grid entites in "
+ << timer.stop() << " seconds." << std::endl;
+ }
+
+ //! the entity set this tree was built with
+ const EntitySet& entitySet() const
+ { return *entitySet_; }
+
+ /////////////////////////////////////////////////////
+ //! Interface to be used by other bounding box trees
+ /////////////////////////////////////////////////////
+
+ //! Get an existing bounding box for a given node
+ const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const
+ { return boundingBoxNodes_[nodeIdx]; }
+
+ //! Get an existing bounding box for a given node
+ const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const
+ { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }
+
+ //! Get the number of bounding boxes currently in the tree
+ std::size_t numBoundingBoxes() const
+ { return boundingBoxNodes_.size(); }
+
+ //! Check whether a bounding box node is a leaf node
+ //! Leaf nodes have itself as child0
+ bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const
+ { return node.child0 == nodeIdx; }
+
+private:
+
+ //! vector of bounding box nodes
+ std::vector<BoundingBoxNode> boundingBoxNodes_;
+
+ //! vector of bounding box coordinates
+ std::vector<ctype> boundingBoxCoordinates_;
+
+ //! a pointer to the entity set
+ std::shared_ptr<const EntitySet> entitySet_;
+
+ //! Compute the bounding box of a grid entity
+ template <class Entity>
+ void computeEntityBoundingBox_(ctype* b, const Entity& entity) const
+ {
+ // get the bounding box coordinates
+ ctype* xMin = b;
+ ctype* xMax = b + dimworld;
+
+ // get mesh entity data
+ auto geometry = entity.geometry();
+
+ // Get coordinates of first vertex
+ auto corner = geometry.corner(0);
+ for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
+ xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
+
+ // Compute the min and max over the remaining vertices
+ for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
+ {
+ corner = geometry.corner(cornerIdx);
+ for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
+ {
+ using std::max;
+ using std::min;
+ xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
+ xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
+ }
+ }
+ }
+
+ //! Build bounding box tree for all entities recursively
+ std::size_t build_(const std::vector<ctype>& leafBoxes,
+ const std::vector<std::size_t>::iterator& begin,
+ const std::vector<std::size_t>::iterator& end)
+ {
+ assert(begin < end);
+
+ // If we reached the end of the recursion, i.e. only a leaf box is left
+ if (end - begin == 1)
+ {
+ // Get the bounding box coordinates for the leaf
+ const std::size_t leafNodeIdx = *begin;
+ const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
+ const auto endCoords = beginCoords + 2*dimworld;
+
+ // Store the data in the bounding box
+ // leaf nodes are indicated by setting child0 to
+ // the node itself and child1 to the index of the entity in the bounding box.
+ return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
+ }
+
+ // Compute the bounding box of all bounding boxes in the range [begin, end]
+ const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
+
+ // sort bounding boxes along the longest axis
+ const auto axis = computeLongestAxis_(bCoords);
+
+ // nth_element sorts the range to make sure that middle points to the coordinate median in axis direction
+ // this is the most expensive part of the algorithm
+ auto middle = begin + (end - begin)/2;
+ std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
+ {
+ const ctype* bi = leafBoxes.data() + 2*dimworld*i;
+ const ctype* bj = leafBoxes.data() + 2*dimworld*j;
+ return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld];
+ });
+
+ // split the bounding boxes into two at the middle iterator and call build recursively, each
+ // call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process.
+ return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
+ bCoords.begin(), bCoords.end());
+ }
+
+ //! Add a new bounding box to the tree
+ template <class Iterator>
+ std::size_t addBoundingBox_(BoundingBoxNode&& node,
+ const Iterator& coordBegin,
+ const Iterator& coordEnd)
+ {
+ // Add the bounding box
+ boundingBoxNodes_.emplace_back(node);
+
+ // Add the bounding box's coordinates
+ boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
+
+ // return the index of the new node
+ return boundingBoxNodes_.size() - 1;
+ }
+
+ //! Compute the bounding box of a vector of bounding boxes
+ std::array<ctype, 2*dimworld>
+ computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes,
+ const std::vector<std::size_t>::iterator& begin,
+ const std::vector<std::size_t>::iterator& end)
+ {
+ std::array<ctype, 2*dimworld> bBoxCoords;
+
+ // copy the iterator and get coordinates of first box
+ auto it = begin;
+ const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
+
+ for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
+ bBoxCoords[coordIdx] = bFirst[coordIdx];
+
+ // Compute min and max with the remaining boxes
+ for (; it != end; ++it)
+ {
+ const auto* b = leafBoxes.data() + 2*dimworld*(*it);
+ for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
+ if (b[coordIdx] < bBoxCoords[coordIdx])
+ bBoxCoords[coordIdx] = b[coordIdx];
+ for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx)
+ if (b[coordIdx] > bBoxCoords[coordIdx])
+ bBoxCoords[coordIdx] = b[coordIdx];
+ }
+
+ return bBoxCoords;
+ }
+
+ //! Compute the bounding box of a vector of bounding boxes
+ std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords)
+ {
+ std::array<ctype, dimworld> axisLength;
+ for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
+ axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
+
+ return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
+ }
+};
+
+/*!
+ * \brief Check whether a point is intersectin a bounding box
+ * \param point The point
+ * \param b Pointer to bounding box coordinates
+ */
+template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0>
+inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
+{
+ static constexpr ctype eps_ = 1.0e-7;
+ const ctype eps0 = eps_*(b[3] - b[0]);
+ const ctype eps1 = eps_*(b[4] - b[1]);
+ const ctype eps2 = eps_*(b[5] - b[2]);
+ return (b[0] - eps0 <= point[0] && point[0] <= b[3] + eps0 &&
+ b[1] - eps1 <= point[1] && point[1] <= b[4] + eps1 &&
+ b[2] - eps2 <= point[2] && point[2] <= b[5] + eps2);
+}
+
+template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0>
+inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
+{
+ static constexpr ctype eps_ = 1.0e-7;
+ const ctype eps0 = eps_*(b[2] - b[0]);
+ const ctype eps1 = eps_*(b[3] - b[1]);
+ return (b[0] - eps0 <= point[0] && point[0] <= b[2] + eps0 &&
+ b[1] - eps1 <= point[1] && point[1] <= b[3] + eps1);
+}
+
+template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0>
+inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
+{
+ static constexpr ctype eps_ = 1.0e-7;
+ const ctype eps0 = eps_*(b[1] - b[0]);
+ return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
+}
+
+/*!
+ * \brief Check whether a bounding box is intersecting another bounding box
+ * \param a Pointer to first bounding box coordinates
+ * \param b Pointer to second bounding box coordinates
+ */
+template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0>
+inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
+{
+ using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
+ static constexpr ctype eps_ = 1.0e-7;
+ const ctype eps0 = eps_*(b[3] - b[0]);
+ const ctype eps1 = eps_*(b[4] - b[1]);
+ const ctype eps2 = eps_*(b[5] - b[2]);
+ return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
+ b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
+ b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
+
+}
+
+template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0>
+inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
+{
+ using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
+ static constexpr ctype eps_ = 1.0e-7;
+ const ctype eps0 = eps_*(b[2] - b[0]);
+ const ctype eps1 = eps_*(b[3] - b[1]);
+ return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
+ b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
+}
+
+template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0>
+inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
+{
+ using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
+ static constexpr ctype eps_ = 1.0e-1;
+ const ctype eps0 = eps_*(b[1] - b[0]);
+ return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/boundingboxtreeintersection.hh b/dumux/common/geometry/boundingboxtreeintersection.hh
new file mode 100644
index 0000000000..d80b906c45
--- /dev/null
+++ b/dumux/common/geometry/boundingboxtreeintersection.hh
@@ -0,0 +1,93 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief A class storing intersections from intersecting two bounding box trees
+ */
+#ifndef DUMUX_BOUNDING_BOX_TREE_INTERSECTION_HH
+#define DUMUX_BOUNDING_BOX_TREE_INTERSECTION_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/promotiontraits.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \brief An intersection object resulting from the intersection of two bounding box tree primitives
+ *
+ * After is has been found that two leaf bounding boxes intersect a primitive test has to be
+ * performed to see if the actual entities inside the bounding box intersect too. The result
+ * if such an intersection is found as an object of this class containing the indices of the
+ * intersecting entities and the corners of the intersection object.
+ */
+template<class EntitySet0, class EntitySet1>
+class BoundingBoxTreeIntersection
+{
+ enum { dimworld = EntitySet0::dimensionworld };
+ using ctype = typename Dune::PromotionTraits<typename EntitySet0::ctype, typename EntitySet1::ctype>::PromotedType;
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+
+public:
+ explicit BoundingBoxTreeIntersection(std::size_t a,
+ std::size_t b,
+ std::vector<GlobalPosition>&& c)
+ : a_(a)
+ , b_(b)
+ , corners_(std::move(c))
+ {
+ static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
+ "Can only store intersections of entity sets with the same world dimension");
+ }
+
+ //! Get the index of the intersecting entity belonging to this grid
+ std::size_t first() const
+ { return a_; }
+
+ //! Get the index of the intersecting entity belonging to the other grid
+ std::size_t second() const
+ { return b_; }
+
+ //! Get the corners of the intersection geometry
+ std::vector<GlobalPosition> corners() const
+ { return corners_; }
+
+ /*!
+ * \brief Check if the corners of this intersection match with the given corners
+ * \note This is useful to check if the intersection geometry of two intersections coincide.
+ */
+ bool cornersMatch(const std::vector<GlobalPosition>& otherCorners) const
+ {
+ if (otherCorners.size() != corners_.size())
+ return false;
+
+ const auto eps = 1.5e-7*(corners_[1] - corners_[0]).two_norm();
+ for (int i = 0; i < corners_.size(); ++i)
+ if ((corners_[i] - otherCorners[i]).two_norm() > eps)
+ return false;
+
+ return true;
+ }
+
+private:
+ std::size_t a_, b_; //!< Indices of the intersection elements
+ std::vector<GlobalPosition> corners_; //!< the corner points of the intersection geometry
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/geometricentityset.hh b/dumux/common/geometry/geometricentityset.hh
new file mode 100644
index 0000000000..e298cb14d0
--- /dev/null
+++ b/dumux/common/geometry/geometricentityset.hh
@@ -0,0 +1,102 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief An interface for a set of geometric entities
+ * \note This can be used e.g. to contruct a bounding box volume hierarchy of a grid
+ * It defines the minimum requirement for such a set
+ */
+#ifndef DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH
+#define DUMUX_GRIDVIEW_GEOMETRIC_ENTITY_SET_HH
+
+#include <memory>
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dumux/common/entitymap.hh>
+
+namespace Dumux {
+
+template <class GridView, int codim = 0>
+class GridViewGeometricEntitySet
+{
+ using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+ using Entity = typename GridView::template Codim<codim>::Entity;
+
+public:
+ GridViewGeometricEntitySet(const GridView& gridView)
+ : gridView_(gridView)
+ , mapper_(gridView, Dune::mcmgLayout(Dune::Codim<codim>()))
+ , entityMap_(std::make_shared<EntityMap<GridView, codim>>(gridView.grid(), mapper_))
+ {}
+
+ GridViewGeometricEntitySet(const GridView& gridView,
+ const Mapper& mapper,
+ std::shared_ptr<const EntityMap<GridView, codim>> entityMap)
+ : gridView_(gridView)
+ , mapper_(mapper)
+ , entityMap_(entityMap)
+ {}
+
+ /*!
+ * \brief The world dimension of the entity set
+ */
+ enum { dimensionworld = GridView::dimensionworld };
+
+ /*!
+ * \brief the coordinate type
+ */
+ using ctype = typename GridView::ctype;
+
+ /*!
+ * \brief the number of entities in this set
+ */
+ decltype(auto) size() const
+ { return gridView_.size(codim); }
+
+ /*!
+ * \brief begin iterator to enable range-based for iteration
+ */
+ decltype(auto) begin() const
+ { return entities(gridView_, Dune::Codim<codim>()).begin(); }
+
+ /*!
+ * \brief end iterator to enable range-based for iteration
+ */
+ decltype(auto) end() const
+ { return entities(gridView_, Dune::Codim<codim>()).end(); }
+
+ /*!
+ * \brief get an entities index
+ */
+ std::size_t index(const Entity& e) const
+ { return mapper_.index(e); }
+
+ /*!
+ * \brief get an entity from an index
+ */
+ Entity entity(std::size_t index) const
+ { return (*entityMap_)[index]; }
+
+private:
+ GridView gridView_;
+ Mapper mapper_;
+ std::shared_ptr<const EntityMap<GridView, codim>> entityMap_;
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometrycollision.hh b/dumux/common/geometry/geometryintersection.hh
similarity index 65%
rename from dumux/common/geometrycollision.hh
rename to dumux/common/geometry/geometryintersection.hh
index d148a578a3..8095a07bb3 100644
--- a/dumux/common/geometrycollision.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -19,8 +19,8 @@
* \brief A class for collision detection of two geometries
* and computation of intersection corners
*/
-#ifndef DUMUX_GEOMETRY_COLLISION_HH
-#define DUMUX_GEOMETRY_COLLISION_HH
+#ifndef DUMUX_GEOMETRY_INTERSECTION_HH
+#define DUMUX_GEOMETRY_INTERSECTION_HH
#include <dune/common/exceptions.hh>
#include <dune/common/promotiontraits.hh>
@@ -41,34 +41,39 @@ template
int dimworld = Geometry1::coorddimension,
int dim1 = Geometry1::mydimension,
int dim2 = Geometry2::mydimension>
-class GeometryCollision
+class GeometryIntersection
{
- using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
- using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
-
public:
+ using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ using IntersectionType = std::vector<GlobalPosition>;
+
//! Determine if the two geometries intersect and compute the intersection corners
- static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector<GlobalPosition>& intersection)
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
{
- static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
- DUNE_THROW(Dune::NotImplemented, "Geometry collision detection with intersection computation not implemented for dimworld = "
+ static_assert(dimworld == Geometry2::coorddimension, "Can only intersect geometries of same coordinate dimension");
+ DUNE_THROW(Dune::NotImplemented, "Geometry intersection detection with intersection computation not implemented for dimworld = "
<< dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
}
};
//! Geometry collision detection with 3d and 1d geometry in 3d space
template <class Geometry1, class Geometry2>
-class GeometryCollision<Geometry1, Geometry2, 3, 3, 1>
+class GeometryIntersection<Geometry1, Geometry2, 3, 3, 1>
{
- static const int dimworld = 3;
- static const int dim1 = 3;
- static const int dim2 = 1;
- static const int dimis = dim2; // the intersection dimension
- using Scalar = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
- using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim1>;
- using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+ enum { dimworld = 3 };
+ enum { dim1 = 3 };
+ enum { dim2 = 1 };
+ enum { dimis = dim2 }; // the intersection dimension
+
+public:
+ using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ using IntersectionType = std::vector<GlobalPosition>;
- static constexpr Scalar eps_ = 1.5e-7; // base epsilon for floating point comparisons
+private:
+ static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
+ using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
public:
/*!
@@ -80,9 +85,9 @@ public:
* \param intersection If the geometries collide intersection holds the corner points of
* the intersection object in global coordinates.
*/
- static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector<GlobalPosition>& intersection)
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
{
- static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension");
+ static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
const auto a = geo2.corner(0);
const auto b = geo2.corner(1);
@@ -91,24 +96,29 @@ public:
// The initial interval is the whole segment
// afterward we start clipping the interval
// by the planes decribed by the facet
- Scalar tfirst = 0;
- Scalar tlast = 1;
+ ctype tfirst = 0.0;
+ ctype tlast = 1.0;
- std::vector<std::vector<int>> facets;
- // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
- switch (geo1.corners())
+ const std::vector<std::vector<int>> facets = [&]()
{
- // todo compute intersection geometries
- case 8: // hexahedron
- facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
- {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
- break;
- case 4: // tetrahedron
- facets = {{1, 2, 0}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
- break;
- default:
- DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type() << ", "<< geo1.corners() << " corners.");
- }
+ std::vector<std::vector<int>> facets;
+ // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
+ switch (geo1.corners())
+ {
+ // todo compute intersection geometries
+ case 8: // hexahedron
+ facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
+ {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
+ break;
+ case 4: // tetrahedron
+ facets = {{1, 2, 0}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
+ break;
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type() << ", "<< geo1.corners() << " corners.");
+ }
+
+ return facets;
+ }();
for (const auto& f : facets)
{
@@ -120,8 +130,8 @@ public:
auto n = crossProduct(v0, v1);
n /= n.two_norm();
- const Scalar denom = n*d;
- const Scalar dist = n*(a-geo1.corner(f[0]));
+ const ctype denom = n*d;
+ const ctype dist = n*(a-geo1.corner(f[0]));
// if denominator is zero the segment in parallel to
// the plane. If the distance is positive there is no intersection
@@ -133,7 +143,7 @@ public:
}
else // not parallel: compute line-plane intersection
{
- const Scalar t = -dist / denom;
+ const ctype t = -dist / denom;
// if entering half space cut tfirst if t is larger
using std::signbit;
if (signbit(denom))
diff --git a/dumux/common/geometry/intersectingentities.hh b/dumux/common/geometry/intersectingentities.hh
new file mode 100644
index 0000000000..0673c61ea5
--- /dev/null
+++ b/dumux/common/geometry/intersectingentities.hh
@@ -0,0 +1,169 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Algorithms that finds which geometric entites intersect
+ */
+#ifndef DUMUX_INTERSECTING_ENTITIES_HH
+#define DUMUX_INTERSECTING_ENTITIES_HH
+
+#include <cmath>
+#include <type_traits>
+#include <dune/common/fvector.hh>
+#include <dumux/common/math.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
+#include <dumux/common/geometry/intersectspointgeometry.hh>
+#include <dumux/common/geometry/geometryintersection.hh>
+#include <dumux/common/geometry/boundingboxtreeintersection.hh>
+
+namespace Dumux
+{
+
+//! Compute all intersections between entities and a point
+template<class EntitySet, class ctype, int dimworld>
+inline std::vector<std::size_t>
+intersectingEntities(const Dune::FieldVector<ctype, dimworld>& point,
+ const BoundingBoxTree<EntitySet>& tree)
+{
+ // Call the recursive find function to find candidates
+ std::vector<std::size_t> entities;
+ intersectingEntities(point, tree, tree.numBoundingBoxes() - 1, entities);
+ return entities;
+}
+
+//! Compute collisions with point for a node of the bounding box tree
+template<class EntitySet, class ctype, int dimworld>
+void intersectingEntities(const Dune::FieldVector<ctype, dimworld>& point,
+ const BoundingBoxTree<EntitySet>& tree,
+ std::size_t node,
+ std::vector<std::size_t>& entities)
+{
+ // Get the bounding box for the current node
+ const auto& bBox = tree.getBoundingBoxNode(node);
+
+ // if the point is not in the bounding box we can stop
+ if (!intersectsPointBoundingBox(point, tree.getBoundingBoxCoordinates(node))) return;
+
+ // now we know it's inside
+ // if the box is a leaf do the primitive test.
+ else if (tree.isLeaf(bBox, node))
+ {
+ const std::size_t entityIdx = bBox.child1;
+ const auto geometry = tree.entitySet().entity(entityIdx).geometry();
+ // if the primitive is positive it intersects the actual geometry, add the entity to the list
+ if (intersectsPointGeometry(point, geometry))
+ entities.push_back(entityIdx);
+ }
+
+ // No leaf. Check both children nodes.
+ else
+ {
+ intersectingEntities(point, tree, bBox.child0, entities);
+ intersectingEntities(point, tree, bBox.child1, entities);
+ }
+}
+
+//! Compute all intersections between two bounding box trees
+template<class EntitySet0, class EntitySet1>
+inline std::vector<BoundingBoxTreeIntersection<EntitySet0, EntitySet1>>
+intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
+ const BoundingBoxTree<EntitySet1>& treeB)
+{
+ // check if the world dimensions match
+ static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
+ "Can only intersect bounding box trees of same world dimension");
+
+ // Create data structure for return type
+ std::vector<BoundingBoxTreeIntersection<EntitySet0, EntitySet1>> intersections;
+
+ // Call the recursive find function to find candidates
+ intersectingEntities(treeA, treeB,
+ treeA.numBoundingBoxes() - 1,
+ treeB.numBoundingBoxes() - 1,
+ intersections);
+
+ return intersections;
+}
+
+//! Compute all intersections between two bounding box tree nodes recursively
+template<class EntitySet0, class EntitySet1>
+void intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
+ const BoundingBoxTree<EntitySet1>& treeB,
+ std::size_t nodeA, std::size_t nodeB,
+ std::vector<BoundingBoxTreeIntersection<EntitySet0, EntitySet1>>& intersections)
+{
+ // Get the bounding box for the current node
+ const auto& bBoxA = treeA.getBoundingBoxNode(nodeA);
+ const auto& bBoxB = treeB.getBoundingBoxNode(nodeB);
+
+ // if the two bounding boxes of the current nodes don't intersect we can stop searching
+ static constexpr int dimworld = EntitySet0::dimensionworld;
+ if (!intersectsBoundingBoxBoundingBox<dimworld>(treeA.getBoundingBoxCoordinates(nodeA),
+ treeB.getBoundingBoxCoordinates(nodeB)))
+ return;
+
+ // Check if we have a leaf in treeA or treeB
+ const bool isLeafA = treeA.isLeaf(bBoxA, nodeA);
+ const bool isLeafB = treeB.isLeaf(bBoxB, nodeB);
+
+ // If both boxes are leaves do the primitive test
+ if (isLeafA && isLeafB)
+ {
+ const auto eIdxA = bBoxA.child1;
+ const auto eIdxB = bBoxB.child1;
+
+ const auto geometryA = treeA.entitySet().entity(eIdxA).geometry();
+ const auto geometryB = treeB.entitySet().entity(eIdxB).geometry();
+
+ using IntersectionAlgorithm = GeometryIntersection< std::decay_t<decltype(geometryA)>,
+ std::decay_t<decltype(geometryB)> >;
+ typename IntersectionAlgorithm::IntersectionType intersection;
+ if (IntersectionAlgorithm::intersection(geometryA, geometryB, intersection))
+ intersections.emplace_back(eIdxA, eIdxB, std::move(intersection));
+ }
+
+ // if we reached the leaf in treeA, just continue in treeB
+ else if (isLeafA)
+ {
+ intersectingEntities(treeA, treeB, nodeA, bBoxB.child0, intersections);
+ intersectingEntities(treeA, treeB, nodeA, bBoxB.child1, intersections);
+ }
+
+ // if we reached the leaf in treeB, just continue in treeA
+ else if (isLeafB)
+ {
+ intersectingEntities(treeA, treeB, bBoxA.child0, nodeB, intersections);
+ intersectingEntities(treeA, treeB, bBoxA.child1, nodeB, intersections);
+ }
+
+ // we know now that both trees didn't reach the leaf yet so
+ // we continue with the larger tree first (bigger node number)
+ else if (nodeA > nodeB)
+ {
+ intersectingEntities(treeA, treeB, bBoxA.child0, nodeB, intersections);
+ intersectingEntities(treeA, treeB, bBoxA.child1, nodeB, intersections);
+ }
+ else
+ {
+ intersectingEntities(treeA, treeB, nodeA, bBoxB.child0, intersections);
+ intersectingEntities(treeA, treeB, nodeA, bBoxB.child1, intersections);
+ }
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/intersectspointgeometry.hh b/dumux/common/geometry/intersectspointgeometry.hh
new file mode 100644
index 0000000000..4b13fd0359
--- /dev/null
+++ b/dumux/common/geometry/intersectspointgeometry.hh
@@ -0,0 +1,111 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Detect if a point intersects a geometry
+ */
+#ifndef DUMUX_INTERSECTS_POINT_GEOMETRY_HH
+#define DUMUX_INTERSECTS_POINT_GEOMETRY_HH
+
+#include <cmath>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+
+#include "intersectspointsimplex.hh"
+
+namespace Dumux {
+
+//! Find out whether a point is inside a three-dimensional geometry
+template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 3), int> = 0>
+bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
+{
+ // get the g type
+ const auto type = g.type();
+
+ // if it's a tetrahedron we can check directly
+ if (type.isTetrahedron())
+ return intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(3));
+
+ // split hexahedrons into five tetrahedrons
+ else if (type.isHexahedron())
+ {
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(3), g.corner(5))) return true;
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(5), g.corner(6), g.corner(4))) return true;
+ if (intersectsPointSimplex(point, g.corner(5), g.corner(3), g.corner(6), g.corner(7))) return true;
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(3), g.corner(2), g.corner(6))) return true;
+ if (intersectsPointSimplex(point, g.corner(5), g.corner(3), g.corner(0), g.corner(6))) return true;
+ return false;
+ }
+
+ // split pyramids into two tetrahedrons
+ else if (type.isPyramid())
+ {
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(4))) return true;
+ if (intersectsPointSimplex(point, g.corner(1), g.corner(3), g.corner(2), g.corner(4))) return true;
+ return false;
+ }
+
+ // split prisms into three tetrahedrons
+ else if (type.isPrism())
+ {
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2), g.corner(4))) return true;
+ if (intersectsPointSimplex(point, g.corner(3), g.corner(0), g.corner(2), g.corner(4))) return true;
+ if (intersectsPointSimplex(point, g.corner(2), g.corner(5), g.corner(3), g.corner(4))) return true;
+ return false;
+ }
+
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "Intersection for point and geometry type "
+ << type << " in " << dimworld << "-dimensional world.");
+}
+
+//! Find out whether a point is inside a two-dimensional geometry
+template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 2), int> = 0>
+bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
+{
+ // get the g type
+ const auto type = g.type();
+
+ // if it's a triangle we can check directly
+ if (type.isTriangle())
+ return intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(2));
+
+ // split quadrilaterals into two triangles
+ else if (type.isQuadrilateral())
+ {
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(1), g.corner(3))) return true;
+ if (intersectsPointSimplex(point, g.corner(0), g.corner(3), g.corner(2))) return true;
+ return false;
+ }
+
+ else
+ DUNE_THROW(Dune::NotImplemented,
+ "Intersection for point and geometry type "
+ << type << " in " << dimworld << "-dimensional world.");
+}
+
+//! Find out whether a point is inside a one-dimensional geometry
+template <class ctype, int dimworld, class Geometry, typename std::enable_if_t<(Geometry::mydimension == 1), int> = 0>
+bool intersectsPointGeometry(const Dune::FieldVector<ctype, dimworld>& point, const Geometry& g)
+{
+ return intersectsPointSimplex(point, g.corner(0), g.corner(1));
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/intersectspointsimplex.hh b/dumux/common/geometry/intersectspointsimplex.hh
new file mode 100644
index 0000000000..601f189057
--- /dev/null
+++ b/dumux/common/geometry/intersectspointsimplex.hh
@@ -0,0 +1,249 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Detect if a point intersects a simplex (including boundary)
+ */
+#ifndef DUMUX_INTERSECTS_POINT_SIMPLEX_HH
+#define DUMUX_INTERSECTS_POINT_SIMPLEX_HH
+
+#include <cmath>
+#include <dune/common/fvector.hh>
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+//! Find out whether a point is inside a tetrahedron (p0, p1, p2, p3) (dimworld is 3)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1,
+ const Dune::FieldVector<ctype, dimworld>& p2,
+ const Dune::FieldVector<ctype, dimworld>& p3)
+{
+ // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
+ // See also "Real-Time Collision Detection" by Christer Ericson.
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // put the tetrahedron points in an array
+ const GlobalPosition *p[4] = {&p0, &p1, &p2, &p3};
+
+ // iterate over all faces
+ for (int i = 0; i < 4; ++i)
+ {
+ // compute all the vectors from vertex (local index 0) to the other points
+ const GlobalPosition v1 = *p[(i + 1)%4] - *p[i];
+ const GlobalPosition v2 = *p[(i + 2)%4] - *p[i];
+ const GlobalPosition v3 = *p[(i + 3)%4] - *p[i];
+ const GlobalPosition v = point - *p[i];
+ // compute the normal to the facet (cross product)
+ GlobalPosition n1 = crossProduct(v1, v2);
+ n1 /= n1.two_norm();
+ // find out on which side of the plane v and v3 are
+ const auto t1 = n1.dot(v);
+ const auto t2 = n1.dot(v3);
+ // If the point is not exactly on the plane the
+ // points have to be on the same side
+ const auto eps = eps_ * v1.two_norm();
+ if ((t1 > eps || t1 < -eps) && std::signbit(t1) != std::signbit(t2))
+ return false;
+ }
+ return true;
+}
+
+//! Find out whether a point is inside a triangle (p0, p1, p2) (dimworld is 3)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1,
+ const Dune::FieldVector<ctype, dimworld>& p2)
+{
+ // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
+ // See also "Real-Time Collision Detection" by Christer Ericson.
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // compute the vectors of the edges and the vector point-p0
+ const GlobalPosition v1 = p0 - p2;
+ const GlobalPosition v2 = p1 - p0;
+ const GlobalPosition v3 = p2 - p1;
+ const GlobalPosition v = point - p0;
+
+ // compute the normal of the triangle
+ const GlobalPosition n = crossProduct(v1, v2);
+
+ // first check if we are in the plane of the triangle
+ // if not we can return early
+ const double t = v.dot(n);
+ using std::abs;
+ if (abs(t) > v1.two_norm()*eps_) // take |v1| as scale
+ return false;
+
+ // compute the normal to the triangle made of point and first edge
+ // the dot product of this normal and the triangle normal has to
+ // be positive because we defined the edges in the right orientation
+ const GlobalPosition n1 = crossProduct(v, v1);
+ const double t1 = n.dot(n1);
+ if (t1 < 0) return false;
+
+ const GlobalPosition n2 = crossProduct(v, v2);
+ const double t2 = n.dot(n2);
+ if (t2 < 0) return false;
+
+ const GlobalPosition n3 = crossProduct(v, v3);
+ const double t3 = n.dot(n3);
+ if (t3 < 0) return false;
+
+ return true;
+}
+
+//! Find out whether a point is inside a triangle (p0, p1, p2) (dimworld is 2)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 2), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1,
+ const Dune::FieldVector<ctype, dimworld>& p2)
+{
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // Use barycentric coordinates
+ const ctype A = 0.5*(-p1[1]*p2[0] + p0[1]*(p2[0] - p1[0])
+ +p1[0]*p2[1] + p0[0]*(p1[1] - p2[1]));
+ const ctype sign = std::copysign(1.0, A);
+ const ctype s = sign*(p0[1]*p2[0] + point[0]*(p2[1]-p0[1])
+ -p0[0]*p2[1] + point[1]*(p0[0]-p2[0]));
+ const ctype t = sign*(p0[0]*p1[1] + point[0]*(p0[1]-p1[1])
+ -p0[1]*p1[0] + point[1]*(p1[0]-p0[0]));
+ const ctype eps = sign*A*eps_;
+
+ return (s > -eps
+ && t > -eps
+ && (s + t) < 2*A*sign + eps);
+}
+
+//! Find out whether a point is inside an interval (p0, p1) (dimworld is 3)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1)
+{
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // compute the vectors between p0 and the other points
+ const GlobalPosition v1 = p1 - p0;
+ const GlobalPosition v2 = point - p0;
+
+ // check if point and p0 are the same
+ const ctype v1norm = v1.two_norm();
+ const ctype v2norm = v2.two_norm();
+ if (v2norm < v1norm*eps_)
+ return true;
+
+ // if not check if p0 and p1 are the same
+ // then we know that point is not in the interval
+ if (v1norm < eps_)
+ return false;
+
+ // if the cross product is zero the points are on a line
+ const GlobalPosition n = crossProduct(v1, v2);
+
+ // early return if the vector length is larger than zero
+ if (n.two_norm() > v1norm*eps_)
+ return false;
+
+ // we know the points are aligned
+ // if the dot product is positive and the length in range
+ // the point is in the interval
+ return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
+}
+
+//! Find out whether a point is inside an interval (p0, p1) (dimworld is 2)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 2), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1)
+{
+ using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // compute the vectors between p0 and the other points
+ const GlobalPosition v1 = p1 - p0;
+ const GlobalPosition v2 = point - p0;
+
+ // check if point and p0 are the same
+ const ctype v1norm = v1.two_norm();
+ const ctype v2norm = v2.two_norm();
+ if (v2norm < v1norm*eps_)
+ return true;
+
+ // if not check if p0 and p1 are the same
+ // then we know that point is not in the interval
+ if (v1norm < eps_)
+ return false;
+
+ // if the cross product is zero the points are on a line
+ const ctype n = crossProduct(v1, v2);
+
+ // early return if the cross product is larger than zero
+ using std::abs;
+ if (abs(n) > v1norm*eps_)
+ return false;
+
+ // we know the points are aligned
+ // if the dot product is positive and the length in range
+ // the point is in the interval
+ return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
+}
+
+//! Find out whether a point is inside an interval (p0, p1) (dimworld is 1)
+template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 1), int> = 0>
+bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
+ const Dune::FieldVector<ctype, dimworld>& p0,
+ const Dune::FieldVector<ctype, dimworld>& p1)
+{
+ static constexpr ctype eps_ = 1.0e-7;
+
+ // sort the interval so interval[1] is the end and interval[0] the start
+ const ctype *interval[2] = {&p0[0], &p1[0]};
+ if (*interval[0] > *interval[1])
+ std::swap(interval[0], interval[1]);
+
+ const ctype v1 = point[0] - *interval[0];
+ const ctype v2 = *interval[1] - *interval[0]; // always positive
+
+ // the coordinates are the same
+ using std::abs;
+ if (abs(v1) < v2*eps_)
+ return true;
+
+ // the point doesn't coincide with p0
+ // so if p0 and p1 are equal it's not inside
+ if (v2 < 1.0e-30)
+ return false;
+
+ // the point is inside if the length is
+ // smaller than the interval length and the
+ // sign of v1 & v2 are the same
+ using std::signbit;
+ return (!signbit(v1) && abs(v1) < v2*(1.0 + eps_));
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 2b429d5df2..bcda706a4d 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -29,7 +29,9 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
-#include <dumux/common/boundingboxtree.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
+#include <dumux/common/geometry/intersectspointgeometry.hh>
+#include <dumux/common/geometry/intersectingentities.hh>
#include <dumux/discretization/methods.hh>
@@ -67,7 +69,7 @@ public:
//! Constructor for sol dependent point sources, when there is no
// value known at the time of initialization
PointSource(GlobalPosition pos)
- : values_(0), pos_(pos), embeddings_(1) {}
+ : values_(0.0), pos_(pos), embeddings_(1) {}
//! Convenience += operator overload modifying only the values
PointSource& operator+= (Scalar s)
@@ -177,7 +179,7 @@ public:
//! Constructor for sol dependent point sources, when there is no
// value known at the time of initialization
IdPointSource(GlobalPosition pos, IdType id)
- : ParentType(pos, PrimaryVariables(0)), id_(id) {}
+ : ParentType(pos, PrimaryVariables(0.0)), id_(id) {}
//! return the sources identifier
IdType id() const
@@ -232,7 +234,7 @@ public:
// value known at the time of initialization
SolDependentPointSource(GlobalPosition pos,
ValueFunction valueFunction)
- : ParentType(pos, PrimaryVariables(0)), valueFunction_(valueFunction) {}
+ : ParentType(pos, PrimaryVariables(0.0)), valueFunction_(valueFunction) {}
//! an update function called before adding the value
// to the local residual in the problem in scvPointSources
@@ -293,7 +295,7 @@ public:
for (auto&& source : sources)
{
// compute in which elements the point source falls
- std::vector<unsigned int> entities = boundingBoxTree.computeEntityCollisions(source.position());
+ const auto entities = intersectingEntities(source.position(), boundingBoxTree);
// split the source values equally among all concerned entities
source.setEmbeddings(entities.size()*source.embeddings());
// loop over all concernes elements
@@ -302,7 +304,7 @@ public:
if(isBox)
{
// check in which subcontrolvolume(s) we are
- const auto element = boundingBoxTree.entity(eIdx);
+ const auto element = boundingBoxTree.entitySet().entity(eIdx);
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
@@ -311,7 +313,7 @@ public:
std::vector<unsigned int> scvIndices;
for (auto&& scv : scvs(fvGeometry))
{
- if (BoundingBoxTreeHelper<dimworld>::pointInGeometry(scv.geometry(), globalPos))
+ if (intersectsPointGeometry(globalPos, scv.geometry()))
scvIndices.push_back(scv.indexInElement());
}
// for all scvs that where tested positiv add the point sources
diff --git a/dumux/discretization/basefvgridgeometry.hh b/dumux/discretization/basefvgridgeometry.hh
index 3787db93d6..8a674d4412 100644
--- a/dumux/discretization/basefvgridgeometry.hh
+++ b/dumux/discretization/basefvgridgeometry.hh
@@ -27,7 +27,9 @@
#include <dune/grid/common/mcmgmapper.hh>
#include <dumux/common/properties.hh>
-#include <dumux/common/boundingboxtree.hh>
+#include <dumux/common/entitymap.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
+#include <dumux/common/geometry/geometricentityset.hh>
namespace Dumux
{
@@ -44,7 +46,9 @@ class BaseFVGridGeometry
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
- using BoundingBoxTree = Dumux::BoundingBoxTree<GridView>;
+ using ElementMap = EntityMap<GridView, 0>;
+ using ElementSet = GridViewGeometricEntitySet<GridView, 0>;
+ using BoundingBoxTree = Dumux::BoundingBoxTree<ElementSet>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -89,8 +93,9 @@ public:
//! Compute the bouding box of the entire domain, for e.g. setting boundary conditions
computeGlobalBoundingBox_();
- //! reset bounding box tree until requested the next time
- boundingBoxTree_.reset(nullptr);
+ //! reset bounding box tree and the element map until requested the next time
+ boundingBoxTree_.release();
+ elementMap_.reset();
}
/*!
@@ -126,23 +131,27 @@ public:
/*!
* \brief Returns the bounding box tree of the grid
*/
- BoundingBoxTree& boundingBoxTree()
+ const BoundingBoxTree& boundingBoxTree() const
{
if(!boundingBoxTree_)
- boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
+ {
+ elementMap(); // make sure the element map is built
+ boundingBoxTree_ = std::make_unique<BoundingBoxTree>
+ ( std::make_shared<ElementSet>(gridView_, elementMapper(), elementMap_) );
+ }
return *boundingBoxTree_;
}
/*!
- * \brief Returns the bounding box tree of the grid
+ * \brief Returns the element index to element map
*/
- const BoundingBoxTree& boundingBoxTree() const
+ const ElementMap& elementMap() const
{
- if(!boundingBoxTree_)
- boundingBoxTree_ = std::make_unique<BoundingBoxTree>(gridView_);
+ if(!elementMap_)
+ elementMap_ = std::make_shared<ElementMap>(gridView_.grid(), elementMapper_);
- return *boundingBoxTree_;
+ return *elementMap_;
}
/*!
@@ -197,6 +206,9 @@ private:
// the bounding box tree of the grid view for effecient element intersections
mutable std::unique_ptr<BoundingBoxTree> boundingBoxTree_;
+ // a map from element index to elements (needed in the bounding box tree and for assembling cell-centered discretization)
+ mutable std::shared_ptr<ElementMap> elementMap_;
+
// the bounding box of the whole domain
GlobalPosition bBoxMin_;
GlobalPosition bBoxMax_;
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 90f1fa01aa..858b4fec5b 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -22,14 +22,12 @@
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
-#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_BASE_HH
-#define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_BASE_HH
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
-#include <dumux/common/elementmap.hh>
-
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/mpfa/connectivitymap.hh>
@@ -81,7 +79,6 @@ public:
//! Per default, we use the secondary IVs at branching points & boundaries
explicit CCMpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
- , elementMap_(gridView)
, secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching)
{ return is.boundary() || isBranching; } )
{}
@@ -89,7 +86,6 @@ public:
//! Constructor with user-defined indicator function for secondary interaction volumes
explicit CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator)
: ParentType(gridView)
- , elementMap_(gridView)
, secondaryIvIndicator_(indicator)
{}
@@ -126,13 +122,13 @@ public:
* \brief Gets an element from a global element index.
*/
Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
+ { return this->elementMap()[eIdx]; }
/*!
* \brief Gets an element from a sub control volume contained in it.
*/
Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
+ { return this->elementMap()[scv.elementIndex()]; }
/*!
* \brief Returns true if primary interaction volumes are used around a given vertex,
@@ -183,7 +179,6 @@ public:
scvs_.resize(numScvs);
scvfs_.reserve(numScvf);
scvfIndicesOfScv_.resize(numScvs);
- elementMap_.resize(numScvs);
// Some methods require to use a second type of interaction volume, e.g.
// around vertices on the boundary or branching points (surface grids)
@@ -206,9 +201,6 @@ public:
// the element geometry
auto elementGeometry = element.geometry();
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
// The local scvf index set
std::vector<IndexType> scvfIndexSet;
scvfIndexSet.reserve(MpfaHelper::getNumLocalScvfs(elementGeometry.type()));
@@ -415,8 +407,7 @@ public:
{ return ivIndexSets_; }
private:
- // mappers
- Dumux::ElementMap<GridView> elementMap_;
+ // connectivity map for efficient assembly
ConnectivityMap connectivityMap_;
// the finite volume grid geometries
@@ -472,7 +463,6 @@ public:
//! Per default, we use the secondary IVs at branching points & boundaries
explicit CCMpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
- , elementMap_(gridView)
, secondaryIvIndicator_([] (const Element& e, const Intersection& is, bool isBranching)
{ return is.boundary() || isBranching; } )
{}
@@ -480,7 +470,6 @@ public:
//! Constructor with user-defined indicator function for secondary interaction volumes
explicit CCMpfaFVGridGeometry(const GridView& gridView, const SecondaryIvIndicator& indicator)
: ParentType(gridView)
- , elementMap_(gridView)
, secondaryIvIndicator_(indicator)
{}
@@ -517,13 +506,13 @@ public:
* \brief Gets an element from a global element index.
*/
Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
+ { return this->elementMap()[eIdx]; }
/*!
* \brief Gets an element from a sub control volume contained in it.
*/
Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
+ { return this->elementMap()[scv.elementIndex()]; }
/*!
* \brief Returns true if primary interaction volumes are used around a given vertex,
@@ -579,7 +568,6 @@ public:
numScvs_ = numDofs();
scvfIndicesOfScv_.resize(numScvs_);
neighborVolVarIndices_.resize(numScvs_);
- elementMap_.resize(numScvs_);
// Some methods require to use a second type of interaction volume, e.g.
// around vertices on the boundary or branching points (surface grids)
@@ -603,9 +591,6 @@ public:
// the element geometry
auto elementGeometry = element.geometry();
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
// the element-wise index sets for finite volume geometry
const auto numLocalFaces = MpfaHelper::getNumLocalScvfs(elementGeometry.type());
std::vector<IndexType> scvfsIndexSet;
@@ -746,8 +731,7 @@ public:
{ return ivIndexSets_; }
private:
- // mappers
- Dumux::ElementMap<GridView> elementMap_;
+ // connectivity map for efficient assembly
ConnectivityMap connectivityMap_;
// containers storing the global data
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index 8add925436..3373a48799 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
@@ -27,9 +27,6 @@
#include <dune/common/version.hh>
-#include <dumux/common/elementmap.hh>
-#include <dumux/common/boundingboxtree.hh>
-
#include <dumux/discretization/basefvgridgeometry.hh>
#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/connectivitymap.hh>
@@ -74,7 +71,6 @@ public:
//! Constructor
CCTpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
- , elementMap_(gridView)
{}
//! the element mapper is the dofMapper
@@ -106,11 +102,11 @@ public:
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
+ { return this->elementMap()[scv.elementIndex()]; }
// Get an element from a global element index
Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
+ { return this->elementMap()[eIdx]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update()
@@ -122,7 +118,6 @@ public:
scvfs_.clear();
scvfIndicesOfScv_.clear();
flipScvfIndices_.clear();
- elementMap_.clear();
// determine size of containers
IndexType numScvs = numDofs();
@@ -131,7 +126,6 @@ public:
numScvf += element.subEntities(1);
// reserve memory
- elementMap_.resize(numScvs);
scvs_.resize(numScvs);
scvfs_.reserve(numScvf);
scvfIndicesOfScv_.resize(numScvs);
@@ -144,9 +138,6 @@ public:
const auto eIdx = this->elementMapper().index(element);
scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
// the element-wise index sets for finite volume geometry
std::vector<IndexType> scvfsIndexSet;
scvfsIndexSet.reserve(element.subEntities(1));
@@ -294,9 +285,8 @@ private:
DUNE_THROW(Dune::InvalidStateException, "No flipped version of this scvf found!");
}
- // mappers
+ // connectivity map for efficient assembly
ConnectivityMap connectivityMap_;
- Dumux::ElementMap<GridView> elementMap_;
// containers storing the global data
std::vector<SubControlVolume> scvs_;
@@ -332,7 +322,6 @@ public:
//! Constructor
CCTpfaFVGridGeometry(const GridView& gridView)
: ParentType(gridView)
- , elementMap_(gridView)
{}
//! the element mapper is the dofMapper
@@ -364,11 +353,11 @@ public:
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
+ { return this->elementMap()[scv.elementIndex()]; }
// Get an element from a global element index
Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
+ { return this->elementMap()[eIdx]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update()
@@ -376,7 +365,6 @@ public:
ParentType::update();
// clear local data
- elementMap_.clear();
scvfIndicesOfScv_.clear();
neighborVolVarIndices_.clear();
@@ -384,7 +372,6 @@ public:
numScvs_ = numDofs();
numScvf_ = 0;
numBoundaryScvf_ = 0;
- elementMap_.resize(numScvs_);
scvfIndicesOfScv_.resize(numScvs_);
neighborVolVarIndices_.resize(numScvs_);
@@ -393,9 +380,6 @@ public:
{
const auto eIdx = this->elementMapper().index(element);
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
// the element-wise index sets for finite volume geometry
auto numLocalFaces = element.subEntities(1);
std::vector<IndexType> scvfsIndexSet;
@@ -489,8 +473,7 @@ private:
IndexType numScvf_;
IndexType numBoundaryScvf_;
- // mappers
- Dumux::ElementMap<GridView> elementMap_;
+ // connectivity map for efficient assembly
ConnectivityMap connectivityMap_;
// vectors that store the global data
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index d18c28c6bd..e1df6246c7 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -26,9 +26,7 @@
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FVGEOMETRY_HH
#include <dumux/common/properties.hh>
-#include <dumux/common/elementmap.hh>
#include <dumux/discretization/basefvgridgeometry.hh>
-#include <dumux/common/boundingboxtree.hh>
#include <dumux/discretization/staggered/connectivitymap.hh>
namespace Dumux
@@ -72,7 +70,8 @@ class StaggeredFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag
public:
//! Constructor
StaggeredFVGridGeometry(const GridView& gridView)
- : ParentType(gridView), elementMap_(gridView), intersectionMapper_(gridView) {}
+ : ParentType(gridView)
+ , intersectionMapper_(gridView) {}
//! The total number of sub control volumes
std::size_t numScv() const
@@ -110,11 +109,11 @@ public:
// Get an element from a sub control volume contained in it
Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
+ { return this->elementMap()[scv.elementIndex()]; }
// Get an element from a global element index
Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
+ { return this->elementMap()[eIdx]; }
//! update all fvElementGeometries (do this again after grid adaption)
void update()
@@ -123,7 +122,6 @@ public:
scvs_.clear();
scvfs_.clear();
scvfIndicesOfScv_.clear();
- elementMap_.clear();
intersectionMapper_.update();
// determine size of containers
@@ -133,7 +131,6 @@ public:
numScvf += element.subEntities(1);
// reserve memory
- elementMap_.resize(numScvs);
scvs_.resize(numScvs);
scvfs_.reserve(numScvf);
scvfIndicesOfScv_.resize(numScvs);
@@ -152,9 +149,6 @@ public:
scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
- // fill the element map with seeds
- elementMap_[eIdx] = element.seed();
-
// the element-wise index sets for finite volume geometry
std::vector<IndexType> scvfsIndexSet;
scvfsIndexSet.reserve(numLocalFaces);
@@ -243,7 +237,6 @@ private:
// mappers
ConnectivityMap connectivityMap_;
- Dumux::ElementMap<GridView> elementMap_;
IntersectionMapper intersectionMapper_;
std::vector<SubControlVolume> scvs_;
diff --git a/test/common/boundingboxtree/CMakeLists.txt b/test/common/boundingboxtree/CMakeLists.txt
index 3ccdfee921..13754bafd4 100644
--- a/test/common/boundingboxtree/CMakeLists.txt
+++ b/test/common/boundingboxtree/CMakeLists.txt
@@ -1,5 +1,15 @@
# build the tests for the bounding box tree
-dune_add_test(SOURCES test_bboxtree.cc)
+dune_add_test(NAME test_bboxtree_dim1
+ SOURCES test_bboxtree.cc
+ COMPILE_DEFINITIONS WORLD_DIMENSION=1)
+
+dune_add_test(NAME test_bboxtree_dim2
+ SOURCES test_bboxtree.cc
+ COMPILE_DEFINITIONS WORLD_DIMENSION=2)
+
+dune_add_test(NAME test_bboxtree_dim3
+ SOURCES test_bboxtree.cc
+ COMPILE_DEFINITIONS WORLD_DIMENSION=3)
# symlink the input file in the build directory
dune_symlink_to_source_files(FILES "network1d.msh")
@@ -8,3 +18,5 @@ dune_symlink_to_source_files(FILES "network1d.msh")
install(FILES
test_bboxtree.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/boundingboxtree)
+
+set(CMAKE_BUILD_TYPE Release)
diff --git a/test/common/boundingboxtree/test_bboxtree.cc b/test/common/boundingboxtree/test_bboxtree.cc
index 0cdccae1a3..9da087f549 100644
--- a/test/common/boundingboxtree/test_bboxtree.cc
+++ b/test/common/boundingboxtree/test_bboxtree.cc
@@ -1,11 +1,13 @@
#include <config.h>
#include <iostream>
+#include <algorithm>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/yaspgrid.hh>
+#include <dune/common/timer.hh>
#if HAVE_DUNE_FOAMGRID
#include <dune/foamgrid/foamgrid.hh>
@@ -18,7 +20,9 @@
#endif
#include <dumux/common/exceptions.hh>
-#include <dumux/common/boundingboxtree.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
+#include <dumux/common/geometry/geometricentityset.hh>
+#include <dumux/common/geometry/intersectingentities.hh>
namespace Dumux {
@@ -29,50 +33,45 @@ class BBoxTreeTests
using Scalar = typename Grid::ctype;
enum { dimWorld = Grid::dimensionworld };
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+ using EntitySet = GridViewGeometricEntitySet<GridView, 0>;
+ using BoundingBoxTree = Dumux::BoundingBoxTree<EntitySet>;
public:
- int construct(const GridView& gv)
+ int build(const GridView& gv)
{
- std::cout << std::endl
- << "Construct bounding box tree from grid view:" << std::endl
- << "*******************************************"
- << std::endl;
-
- // construct a bounding box tree
- tree_ = std::make_shared<Dumux::BoundingBoxTree<GridView>>();
- tree_->build(gv);
-
+ // build a bounding box tree
+ tree_ = std::make_shared<BoundingBoxTree>();
+ tree_->build(std::make_shared<EntitySet>(gv));
return 0;
}
int intersectPoint(const GlobalPosition& p, std::size_t expectedCollisions)
{
- std::cout << std::endl
- << "Intersect with point ("<< p <<"):" << std::endl
- << "************************************************"
- << std::endl;
+ std::cout << "Intersect with point ("<< p <<") ";
+
+ Dune::Timer timer;
+ const auto entities = intersectingEntities(p, *tree_);
- auto entities = tree_->computeEntityCollisions(p);
+ std::cout << " --> " << entities.size() << " intersection(s) found ("
+ << expectedCollisions << " expected) in " << timer.elapsed() << " seconds.\n";
- std::cout << entities.size() << " intersection(s) found ("
- << expectedCollisions << " expected)." << std::endl;
if (entities.size() != expectedCollisions)
{
std::cerr << "Point intersection failed: Expected "
<< expectedCollisions << " and got "
- << entities.size() << "!" <<std::endl;
+ << entities.size() << "!\n";
return 1;
}
return 0;
}
- template <class OtherGridView>
- int intersectTree(const Dumux::BoundingBoxTree<OtherGridView>& otherTree,
+ template <class OtherEntitySet, class OtherGridView>
+ int intersectTree(const Dumux::BoundingBoxTree<OtherEntitySet>& otherTree,
const OtherGridView& otherGridView,
std::size_t expectedIntersections)
{
Dune::Timer timer;
- auto intersections = tree_->computeEntityCollisions(otherTree);
+ const auto intersections = intersectingEntities(*tree_, otherTree);
std::cout << "Computed tree intersections in " << timer.elapsed() << std::endl;
timer.reset();
@@ -108,14 +107,7 @@ public:
}
private:
- template<class Intersection>
- bool intersectionsEqual(const Intersection& is1, const Intersection& is2)
- {
- using std::max;
- auto eps = 1e-7*max((is1[0] - is1[1]).two_norm(), (is2[0] - is2[1]).two_norm());
- return (is1[0] - is2[0]).two_norm() < eps && (is1[1] - is2[1]).two_norm() < eps;
- }
- std::shared_ptr<Dumux::BoundingBoxTree<GridView>> tree_;
+ std::shared_ptr<BoundingBoxTree> tree_;
};
} // end namespace Dumux
@@ -126,7 +118,8 @@ int main (int argc, char *argv[]) try
Dune::MPIHelper::instance(argc, argv);
// Some aliases two type tags for tests using two grids
- using Grid = Dune::YaspGrid<3>;
+ constexpr int dimworld = WORLD_DIMENSION;
+ using Grid = Dune::YaspGrid<dimworld>;
using Scalar = Grid::ctype;
enum { dimWorld = Grid::dimensionworld };
enum { dim = Grid::dimension };
@@ -144,48 +137,69 @@ int main (int argc, char *argv[]) try
<< std::endl;
// create a cube grid
- const GlobalPosition lowerLeft(0.0);
- const GlobalPosition upperRight(1.0*scaling);
- std::array<unsigned int, dim> elems; elems.fill(3);
- auto grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elems);
+ {
+ const GlobalPosition lowerLeft(0.0);
+ const GlobalPosition upperRight(1.0*scaling);
+ constexpr int numCellsX = 33;
+ std::array<unsigned int, dim> elems; elems.fill(numCellsX);
+ auto grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elems);
+
+ // Dune::VTKWriter<Grid::LeafGridView> vtkWriter(grid->leafGridView());
+ // vtkWriter.write("grid_dim" + std::to_string(dimworld), Dune::VTK::ascii);
+
+ // bboxtree tests using one bboxtree
+ returns.push_back(test.build(grid->leafGridView()));
+ returns.push_back(test.intersectPoint(GlobalPosition(0.0*scaling), 1));
+ returns.push_back(test.intersectPoint(GlobalPosition(1e-3*scaling), 1));
+ returns.push_back(test.intersectPoint(GlobalPosition(1.0*scaling/numCellsX), 1<<dimworld));
+ returns.push_back(test.intersectPoint(GlobalPosition(1.0*scaling), 1));
+ }
+
+#if HAVE_DUNE_FOAMGRID && WORLD_DIMENSION == 3
+ {
+ const GlobalPosition lowerLeft(0.0);
+ const GlobalPosition upperRight(1.0*scaling);
+ constexpr int numCellsX = 10;
+ std::array<unsigned int, dim> elems; elems.fill(numCellsX);
+ auto grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lowerLeft, upperRight, elems);
- Dune::VTKWriter<Grid::LeafGridView> vtkWriter(grid->leafGridView());
- vtkWriter.write("grid", Dune::VTK::ascii);
+ // Dune::VTKWriter<Grid::LeafGridView> vtkWriter(grid->leafGridView());
+ // vtkWriter.write("grid_dim" + std::to_string(dimworld), Dune::VTK::ascii);
- // bboxtree tests using one bboxtree
- returns.push_back(test.construct(grid->leafGridView()));
- returns.push_back(test.intersectPoint(GlobalPosition(0.0), 1));
- returns.push_back(test.intersectPoint(GlobalPosition(1e-3*scaling), 1));
- returns.push_back(test.intersectPoint(GlobalPosition(1.0/3.0*scaling), 8));
+ using NetworkGrid = Dune::FoamGrid<1, dimworld>;
+ using NetworkGridView = NetworkGrid::LeafGridView;
-#if HAVE_DUNE_FOAMGRID
- using NetworkGrid = Dune::FoamGrid<1, 3>;
- using NetworkGridView = NetworkGrid::LeafGridView;
+ std::cout << std::endl
+ << "Intersect with other bounding box tree:" << std::endl
+ << "***************************************"
+ << std::endl;
- std::cout << std::endl
- << "Intersect with other bounding box tree:" << std::endl
- << "***************************************"
- << std::endl;
+ auto networkGrid = std::shared_ptr<NetworkGrid>(Dune::GmshReader<NetworkGrid>::read("network1d.msh", false, false));
- auto networkGrid = std::shared_ptr<NetworkGrid>(Dune::GmshReader<NetworkGrid>::read("network1d.msh", false, false));
+ // scaling
+ for (const auto& vertex : vertices(networkGrid->leafGridView()))
+ {
+ auto newPos = vertex.geometry().corner(0);
+ newPos *= scaling;
+ networkGrid->setPosition(vertex, newPos);
+ }
- // scaling
- for (const auto& vertex : vertices(networkGrid->leafGridView()))
- {
- auto newPos = vertex.geometry().corner(0);
- newPos *= scaling;
- networkGrid->setPosition(vertex, newPos);
- }
+ // Dune::VTKWriter<NetworkGridView> lowDimVtkWriter(networkGrid->leafGridView());
+ // lowDimVtkWriter.write("network", Dune::VTK::ascii);
- Dune::VTKWriter<NetworkGridView> lowDimVtkWriter(networkGrid->leafGridView());
- lowDimVtkWriter.write("network", Dune::VTK::ascii);
+ std::cout << "Constructed 1d network grid with " << networkGrid->leafGridView().size(0) << " elements." << std::endl;
- std::cout << "Constructed " << networkGrid->leafGridView().size(0) <<
- " element 1d network grid." << std::endl;
+ // build the bulk grid bounding box tree
+ returns.push_back(test.build(grid->leafGridView()));
- Dumux::BoundingBoxTree<NetworkGridView> networkTree;
- networkTree.build(networkGrid->leafGridView());
- returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 9));
+ // build the network grid bounding box tree
+ using EntitySet = Dumux::GridViewGeometricEntitySet<NetworkGridView, 0>;
+ Dumux::BoundingBoxTree<EntitySet> networkTree;
+ networkTree.build(std::make_shared<EntitySet>(networkGrid->leafGridView()));
+
+ // intersect the two bounding box trees
+ returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 20));
+ }
#endif
}
--
GitLab