diff --git a/dumux/common/boundingboxtree.hh b/dumux/common/boundingboxtree.hh index e00a21acbf1ffcbf864137985a5fecfaeadbb6d7..50de581baf63007c9dc152a24814403dd24f6236 100644 --- a/dumux/common/boundingboxtree.hh +++ b/dumux/common/boundingboxtree.hh @@ -16,9 +16,9 @@ *****************************************************************************/ /*! * \file - * \brief A bounding box tree for grid point intersections + * \brief An axis-aligned bounding box volume hierarchy for dune grids * - * Dumux implementation of the bounding box tree + * Dumux implementation of an AABB tree * adapted from implementation in FEniCS by Anders Logg */ #ifndef DUMUX_BOUNDINGBOXTREE_HH @@ -30,13 +30,17 @@ #include #include #include +#include namespace Dumux { -// optimized dimension-dependent methods +/*! + * \brief A helper class to optimize dimension-dependent methods + */ template class BoundingBoxTreeHelper {}; +//! Helper methods for world dimension three template <> class BoundingBoxTreeHelper<3> { @@ -45,7 +49,7 @@ class BoundingBoxTreeHelper<3> typedef Dune::FieldVector GlobalPosition; public: - // Check whether a point is inside a given geometry + //! Check whether a point is inside a given three-dimensional geometry template static typename std::enable_if::type pointInGeometry(const Geometry& geometry, const GlobalPosition& point) @@ -143,6 +147,7 @@ public: << type << " in three-dimensional world."); } + //! Check whether a point is inside a given two-dimensional geometry template static typename std::enable_if::type pointInGeometry(const Geometry& geometry, const GlobalPosition& point) @@ -180,6 +185,7 @@ public: } + //! Check whether a point is inside a given one-dimensional geometry template static typename std::enable_if::type pointInGeometry(const Geometry& geometry, const GlobalPosition& point) @@ -189,7 +195,7 @@ public: point); } - //! find out if a point is inside a tetrahedron + //! 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) @@ -210,20 +216,21 @@ public: const GlobalPosition v3 = *p[(i + 3)%4] - *p[i]; const GlobalPosition v = point - *p[i]; // compute the normal to the facet (cross product) - const GlobalPosition n1 = Dumux::crossProduct(v1, v2); + GlobalPosition n1 = Dumux::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_ * v.two_norm(); + const double eps = eps_ * v1.two_norm(); if ((t1 > eps || t1 < -eps) && std::signbit(t1) != std::signbit(t2)) return false; } return true; } - //! find out if a point is inside a triangle + //! 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) { @@ -264,7 +271,7 @@ public: return true; } - //! find out if a point is inside an interval + //! Find out whether a point is inside an interval (p0, p1) static bool pointInInterval(const GlobalPosition& p0, const GlobalPosition& p1, const GlobalPosition& point) { @@ -299,12 +306,14 @@ public: return false; } - // Check whether a point is in a bounding box + /*! + * \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& point, - const double* boundingBoxCoordinates, - unsigned int node) + const double* b) { - const double* b = boundingBoxCoordinates + 6*node; const double eps0 = eps_*(b[3] - b[0]); const double eps1 = eps_*(b[4] - b[1]); const double eps2 = eps_*(b[5] - b[2]); @@ -313,7 +322,22 @@ public: b[2] - eps2 <= point[2] && point[2] <= b[5] + eps2); } - // Compute the bounding box of a vector of bounding boxes + /*! + * \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& leafBoxes, @@ -352,7 +376,7 @@ public: axis = 2; } - // Sort the bounding boxes along the longest axis + //! Sort the bounding boxes along the longest axis static void sortBoundingBoxes(std::size_t axis, const std::vector& leafBoxes, const std::vector::iterator& begin, @@ -370,8 +394,10 @@ public: } } - // Comparison operators for sorting bounding boxes. There are sorted by their - // mid points along the longest axis + /*! + * \brief Comparison function for sorting bounding boxes on the x-axis + * \note This could be replaced by lambdas + */ struct lessXBox { const std::vector& bBoxes; @@ -384,6 +410,10 @@ public: } }; + /*! + * \brief Comparison function for sorting bounding boxes on the y-axis + * \note This could be replaced by lambdas + */ struct lessYBox { const std::vector& bBoxes; @@ -396,6 +426,10 @@ public: } }; + /*! + * \brief Comparison function for sorting bounding boxes on the z-axis + * \note This could be replaced by lambdas + */ struct lessZBox { const std::vector& bBoxes; @@ -409,6 +443,7 @@ public: }; }; +//! Helper methods for world dimension two template <> class BoundingBoxTreeHelper<2> { @@ -517,19 +552,34 @@ public: return false; } - // Check whether a point is in a bounding box + /*! + * \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& point, - const double* boundingBoxCoordinates, - unsigned int node) + const double* b) { - const double* b = boundingBoxCoordinates + 4*node; 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); } - // Compute the bounding box of a vector of bounding boxes + /*! + * \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& leafBoxes, @@ -563,7 +613,7 @@ public: axis = 1; } - // Sort the bounding boxes along the longest axis + //! Sort the bounding boxes along the longest axis static void sortBoundingBoxes(std::size_t axis, const std::vector& leafBoxes, const std::vector::iterator& begin, @@ -576,8 +626,10 @@ public: std::nth_element(begin, middle, end, lessYBox(leafBoxes)); } - // Comparison operators for sorting bounding boxes. There are sorted by their - // mid points along the longest axis + /*! + * \brief Comparison function for sorting bounding boxes on the x-axis + * \note This could be replaced by lambdas + */ struct lessXBox { const std::vector& bBoxes; @@ -590,6 +642,10 @@ public: } }; + /*! + * \brief Comparison function for sorting bounding boxes on the y-axis + * \note This could be replaced by lambdas + */ struct lessYBox { const std::vector& bBoxes; @@ -603,6 +659,7 @@ public: }; }; +//! Helper methods for world dimension one template <> class BoundingBoxTreeHelper<1> { @@ -644,17 +701,30 @@ public: return false; } - // Check whether a point is in a bounding box + /*! + * \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* boundingBoxCoordinates, - unsigned int node) + const double* b) { - const double* b = boundingBoxCoordinates + 2*node; const double eps0 = eps_*(b[1] - b[0]); return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0; } - // Compute the bounding box of a vector of bounding boxes + /*! + * \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& leafBoxes, @@ -679,7 +749,7 @@ public: axis = 0; } - // Sort the bounding boxes along the longest axis + //! Sort the bounding boxes along the longest axis static void sortBoundingBoxes(std::size_t axis, const std::vector& leafBoxes, const std::vector::iterator& begin, @@ -687,8 +757,10 @@ public: const std::vector::iterator& end) { std::nth_element(begin, middle, end, lessXBox(leafBoxes)); } - // Comparison operators for sorting bounding boxes. There are sorted by their - // mid points along the longest axis + /*! + * \brief Comparison function for sorting bounding boxes on the x-axis + * \note This could be replaced by lambdas + */ struct lessXBox { const std::vector& bBoxes; @@ -702,7 +774,62 @@ public: }; }; -//! An index to element map +/*! + * \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 BoundingBoxTreeIntersection +{ + const static int dimworld = GridView1::dimensionworld; + using Scalar = typename GridView1::ctype; + using GlobalPosition = Dune::FieldVector; + +public: + BoundingBoxTreeIntersection(unsigned int a, unsigned int b, std::vector&& 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 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& 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 corners_; //!< the corner points of the intersection geometry +}; + +/*! + * \brief A class mapping from an element index to elements using element seeds + */ template class IndexToElementMap : public std::vector::EntitySeed> @@ -713,6 +840,7 @@ public: IndexToElementMap(const GridView& gridView) : grid_(gridView.grid()) {} + //! get an element from an index t template Element entity(T&& t) { return grid_.entity((*this)[std::forward(t)]); } @@ -721,28 +849,36 @@ private: const Grid& grid_; }; -//! The bounding box class. Implements an axis-aligned bounding box tree for grids. +/*! + * \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 BoundingBoxTree { + //! be friends with all other kinds bounding box trees so that + //! they can call each others private methods + template friend class BoundingBoxTree; + static const int dim = GridView::dimension; static const int dimworld = GridView::dimensionworld; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename Dune::ReferenceElements ReferenceElements; - typedef typename Dune::ReferenceElement ReferenceElement; + using Element = typename GridView::template Codim<0>::Entity; + using Scalar = typename GridView::ctype; + using GlobalPosition = typename Dune::FieldVector; public: - // Default Constructor + //! Default Constructor BoundingBoxTree() {} - // Constructor with gridView + //! Constructor with gridView BoundingBoxTree(const GridView& leafGridView) { build(leafGridView); } - // Destructor - ~BoundingBoxTree() {} - - // Build up bounding box tree for a grid + //! Build up bounding box tree for a grid with leafGridView void build(const GridView& leafGridView) { // clear data if any @@ -782,7 +918,7 @@ public: << timer.stop() << " seconds." << std::endl; } - // Compute all intersections between entities and a point + //! Compute all intersections between entities and a point std::vector computeEntityCollisions(const Dune::FieldVector& point) const { // Call the recursive find function to find candidates @@ -791,30 +927,52 @@ public: return entities; } - // Get an element from a global element index + //! Compute all intersections between entities and another bounding box tree + template + std::vector> + computeEntityCollisions(const Dumux::BoundingBoxTree& 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> 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: - - // Bounding box data. 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. + /*! + * \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 + //! Vector of bounding boxes std::vector boundingBoxes_; - // Vector of bounding box coordinates + //! Vector of bounding box coordinates std::vector boundingBoxCoordinates_; - // Shared pointer to the index to element map + //! Shared pointer to the index to element map std::shared_ptr > indexToElementMap_; - // Clear all data + //! Clear all data void clear_() { boundingBoxes_.clear(); @@ -822,7 +980,7 @@ private: if(indexToElementMap_) indexToElementMap_->clear(); } - // Build bounding box tree for all entities recursively + //! Build bounding box tree for all entities recursively unsigned int build_(const std::vector& leafBoxes, const std::vector::iterator& begin, const std::vector::iterator& end) @@ -862,7 +1020,7 @@ private: return addBoundingBox_(bBox, b); } - // Compute collisions with point recursively + //! Compute collisions with point recursively void computeCollisions_(const Dune::FieldVector& point, unsigned int node, std::vector& entities) const @@ -871,7 +1029,7 @@ private: const BoundingBox& bBox = getBoundingBox_(node); // if the point is not in the bounding box we can stop - if (!BoundingBoxTreeHelper::pointInBoundingBox(point, boundingBoxCoordinates_.data(), node)) + if (!BoundingBoxTreeHelper::pointInBoundingBox(point, getBoundingBoxCoordinates_(node))) return; // We know now it's inside. If the box is a leaf add it. @@ -896,7 +1054,74 @@ private: } } - // Add a new bounding box to the tree + //! Compute collisions with other bounding box tree recursively + template + void computeCollisions_(const Dumux::BoundingBoxTree& treeB, + unsigned int nodeA, + unsigned int nodeB, + std::vector>& 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:: + 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 = Dumux::GeometryCollision; + std::vector 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) { @@ -911,20 +1136,24 @@ private: return boundingBoxes_.size() - 1; } - // Get an existing bounding box for a given node + //! Get an existing bounding box for a given node inline const BoundingBox& getBoundingBox_(unsigned int node) const { return boundingBoxes_[node]; } - // Get the number of bounding boxes currently in the tree + //! 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 + //! 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 + //! Compute the bounding box of a grid entity template void computeEntityBoundingBox_(double* b, const Entity& entity) const { diff --git a/dumux/common/geometrycollision.hh b/dumux/common/geometrycollision.hh new file mode 100644 index 0000000000000000000000000000000000000000..a7d5e1823e8252d4dcce5f2dd3ed41f2d88a33b7 --- /dev/null +++ b/dumux/common/geometrycollision.hh @@ -0,0 +1,164 @@ +/***************************************************************************** + * 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 . * + *****************************************************************************/ +/*! + * \file + * \brief A class for collision detection of two geometries + * and computation of intersection corners + */ +#ifndef DUMUX_GEOMETRY_COLLISION_HH +#define DUMUX_GEOMETRY_COLLISION_HH + +#include +#include +#include + +#include + +namespace Dumux +{ + +/*! + * \brief A class for geometry collision detection and intersection calculation + * The class can be specialized for combinations of dimworld, dim1, dim2, where + * dimworld is the world dimension embedding a grid of dim1 and a grid of dim2. + */ +template + +class GeometryCollision +{ + using Scalar = typename Dune::PromotionTraits::PromotedType; + using GlobalPosition = Dune::FieldVector; + +public: + //! Determine if the two geometries intersect and compute the intersection corners + static bool collide(const Geometry1& geo1, const Geometry2& geo2, std::vector& 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 = " + << dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2); + } +}; + +//! Geometry collision detection with 3d and 1d geometry in 3d space +template +class GeometryCollision +{ + 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::PromotedType; + using ReferenceElements = typename Dune::ReferenceElements; + using GlobalPosition = Dune::FieldVector; + + static constexpr Scalar eps_ = 1.5e-7; // base epsilon for floating point comparisons + +public: + /*! + * \brief Colliding segment and convex polyhedron + * \note Algorithm from "Real-Time Collision Detection" by Christer Ericson + * Basis is the theorem that for any two non-intersecting convex polyhedrons + * a separating plane exists. + * \param geo1/geo2 The geometries to intersect + * \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& intersection) + { + static_assert(dimworld == Geometry2::coorddimension, "Can only collide geometries of same coordinate dimension"); + + const auto a = geo2.corner(0); + const auto b = geo2.corner(1); + const auto d = b - a; + + // 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; + + std::vector> 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."); + } + + for (const auto& f : facets) + { + // compute normal vector by cross product + const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]); + const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]); + const auto eps = eps_*v0.two_norm(); + + auto n = Dumux::crossProduct(v0, v1); + n /= n.two_norm(); + + const Scalar denom = n*d; + const Scalar 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 + if (std::abs(denom) < eps) + { + if (dist > eps) + return false; + } + else // not parallel: compute line-plane intersection + { + const Scalar t = -dist / denom; + // if entering half space cut tfirst if t is larger + if (std::signbit(denom)) + { + if (t > tfirst) + tfirst = t; + } + // if exiting half space cut tlast if t is smaller + else + { + if (t < tlast) + tlast = t; + } + // there is no intersection of the interval is empty + // use unscaled epsilon since t is in local coordinates + if (tfirst > tlast - eps_) + return false; + + } + } + // If we made it until here an intersection exists. We also export + // the intersections geometry now s(t) = a + t(b-a) in [tfirst, tlast] + intersection = {geo2.global(tfirst), geo2.global(tlast)}; + return true; + } +}; + +} // end namespace Dumux + +# endif diff --git a/test/common/CMakeLists.txt b/test/common/CMakeLists.txt index e1d3eb19eb9af2d129a3e98029bfcc4e3d346431..1fed13a61f8032b1e507dbdebd35ad021245ff5e 100644 --- a/test/common/CMakeLists.txt +++ b/test/common/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory("generalproblem") add_subdirectory("propertysystem") add_subdirectory("spline") +add_subdirectory("boundingboxtree") diff --git a/test/common/boundingboxtree/CMakeLists.txt b/test/common/boundingboxtree/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..d42573b7c49176771aaab8b929569412b90197b2 --- /dev/null +++ b/test/common/boundingboxtree/CMakeLists.txt @@ -0,0 +1,12 @@ +# build the tests for the bounding box tree +add_dumux_test(test_bboxtree test_bboxtree test_bboxtree.cc + ${CMAKE_CURRENT_BINARY_DIR}/test_bboxtree) + +# symlink the input file in the build directory +dune_symlink_to_source_files(FILES "network1d.msh") + +#install sources +install(FILES +test_bboxtree.cc +test_bboxtree.input +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/boundingboxtree) \ No newline at end of file diff --git a/test/common/boundingboxtree/network1d.geo b/test/common/boundingboxtree/network1d.geo new file mode 100644 index 0000000000000000000000000000000000000000..300a103fca8eb8fae16e25998ff7365f07a6c3fb --- /dev/null +++ b/test/common/boundingboxtree/network1d.geo @@ -0,0 +1,17 @@ +cl_ = 0.5; + +Point(1) = {0.5, 0.5, 0.5, cl_}; +Point(2) = {0.5, 0.5, 0.9, cl_}; +Point(3) = {0.1, 0.1, 0.1, cl_}; +Point(4) = {0.1, 0.9, 0.1, cl_}; + +Point(5) = {0.9, 0.9, 0.1, cl_}; +Point(6) = {0.9, 0.9, 0.9, cl_}; + +Line(1) = {2, 1}; +Line(2) = {1, 3}; +Line(3) = {1, 4}; +Line(4) = {5, 6}; + +Transfinite Line{1, 2, 3} = 1; +Transfinite Line{4} = 2; diff --git a/test/common/boundingboxtree/network1d.msh b/test/common/boundingboxtree/network1d.msh new file mode 100644 index 0000000000000000000000000000000000000000..87274e3ca40663592537664eba4a78038e16867c --- /dev/null +++ b/test/common/boundingboxtree/network1d.msh @@ -0,0 +1,25 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +6 +1 0.5 0.5 0.5 +2 0.5 0.5 0.9 +3 0.1 0.1 0.1 +4 0.1 0.9 0.1 +5 0.9 0.9 0.1 +6 0.9 0.9 0.9 +$EndNodes +$Elements +10 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 1 2 0 1 2 1 +8 1 2 0 2 1 3 +9 1 2 0 3 1 4 +10 1 2 0 4 5 6 +$EndElements diff --git a/test/common/boundingboxtree/test_bboxtree.cc b/test/common/boundingboxtree/test_bboxtree.cc new file mode 100644 index 0000000000000000000000000000000000000000..fa5cb0870bb98165e18cf942a6eed8fe5d2b9f35 --- /dev/null +++ b/test/common/boundingboxtree/test_bboxtree.cc @@ -0,0 +1,212 @@ +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace Dumux { + +namespace Properties +{ + NEW_TYPE_TAG(BBoxTreeTest, INHERITS_FROM(NumericModel)); +#if HAVE_UG + SET_TYPE_PROP(BBoxTreeTest, Grid, Dune::UGGrid<3>); +#else + SET_TYPE_PROP(BBoxTreeTest, Grid, Dune::YaspGrid<3>); +#endif +} + +template +class BBoxTreeTests +{ + using Grid = typename GET_PROP_TYPE(TypeTag, Grid); + using GridView = typename Grid::LeafGridView; + using Scalar = typename Grid::ctype; + enum { dimWorld = Grid::dimensionworld }; + using GlobalPosition = Dune::FieldVector; + +public: + int construct(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>(); + tree_->build(gv); + + return 0; + } + + int intersectPoint(const GlobalPosition& p, std::size_t expectedCollisions) + { + std::cout << std::endl + << "Intersect with point ("<< p <<"):" << std::endl + << "************************************************" + << std::endl; + + auto entities = tree_->computeEntityCollisions(p); + + 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() << "!" < + int intersectTree(const Dumux::BoundingBoxTree& otherTree, + const OtherGridView& otherGridView, + std::size_t expectedIntersections) + { + Dune::Timer timer; + auto intersections = tree_->computeEntityCollisions(otherTree); + std::cout << "Computed tree intersections in " << timer.elapsed() << std::endl; + timer.reset(); + + std::vector>> map; + map.resize(otherGridView.size(0)); + std::vector> uniqueIntersections; + for (const auto& is : intersections) + { + bool add = true; + for (const auto& i : map[is.second()]) + { + add = !is.cornersMatch(i); + if (!add) break; + } + if(add) + { + map[is.second()].push_back(is.corners()); + uniqueIntersections.push_back(is.corners()); + } + } + + std::cout << "Found " << uniqueIntersections.size() << " unique intersections " + << "in " << timer.elapsed() << std::endl; + + if (uniqueIntersections.size() != expectedIntersections) + { + std::cerr << "BoundingBoxTree intersection failed: Expected " + << expectedIntersections << " and got " + << uniqueIntersections.size() << "!" < + bool intersectionsEqual(const Intersection& is1, const Intersection& is2) + { + auto eps = 1e-7*std::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> tree_; +}; + +} // end namespace Dumux + +int main (int argc, char *argv[]) try +{ + // maybe initialize mpi + Dune::MPIHelper::instance(argc, argv); + + // Some aliases two type tags for tests using two grids + using TypeTag = TTAG(BBoxTreeTest); + using Grid = GET_PROP_TYPE(TypeTag, Grid); + using Scalar = Grid::ctype; + enum { dimWorld = Grid::dimensionworld }; + enum { dim = Grid::dimension }; + using GlobalPosition = Dune::FieldVector; + + // collect returns to determine exit code + std::vector returns; + Dumux::BBoxTreeTests test; + + for (const auto scaling : {1e10, 1.0, 1e-3, 1e-10}) + { + std::cout << std::endl + << "Testing with scaling = " << scaling << std::endl + << "***************************************" + << std::endl; + + // create a cube grid + const GlobalPosition lowerLeft(0.0); + const GlobalPosition upperRight(1.0*scaling); + std::array elems; elems.fill(3); + auto grid = Dune::StructuredGridFactory::createCubeGrid(lowerLeft, upperRight, elems); + + Dune::VTKWriter vtkWriter(grid->leafGridView()); + vtkWriter.write("grid", 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)); + +#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; + + auto networkGrid = std::shared_ptr(Dune::GmshReader::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); + } + + Dune::VTKWriter lowDimVtkWriter(networkGrid->leafGridView()); + lowDimVtkWriter.write("network", Dune::VTK::ascii); + + std::cout << "Constructed " << networkGrid->leafGridView().size(0) << + " element 1d network grid." << std::endl; + + Dumux::BoundingBoxTree networkTree; + networkTree.build(networkGrid->leafGridView()); + returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 9)); +#endif + } + + std::cout << std::endl; + + // determine the exit code + if (std::any_of(returns.begin(), returns.end(), [](int i){ return i==1; })) + return 1; + + return 0; +} +// ////////////////////////////////// +// Error handler +// ///////////////////////////////// +catch (Dumux::ParameterException &e) { + std::cerr << e << ". Abort!\n"; + return 1; +} +catch (const Dune::Exception& e) { + std::cout << e << std::endl; + return 1; +}