diff --git a/dumux/common/geometry/boundingboxtree.hh b/dumux/common/geometry/boundingboxtree.hh
index 3541867287b4d2425bafdce860e2e313daf5d201..9ca148086ac58187cbb96f0f75dca404da2cc74c 100644
--- a/dumux/common/geometry/boundingboxtree.hh
+++ b/dumux/common/geometry/boundingboxtree.hh
@@ -18,399 +18,15 @@
* \file
* \ingroup Geometry
* \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.
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_BOUNDINGBOXTREE_HH
-#define DUMUX_BOUNDINGBOXTREE_HH
+#ifndef DUMUX_COMMON_GEOMETRY_BOUNDINGBOXTREE_HH
+#define DUMUX_COMMON_GEOMETRY_BOUNDINGBOXTREE_HH
-#include <vector>
-#include <array>
-#include <algorithm>
-#include <memory>
-#include <numeric>
-#include <type_traits>
-#include <iostream>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/boundingboxtree.hh"
-#include <dune/common/promotiontraits.hh>
-#include <dune/common/timer.hh>
-#include <dune/common/fvector.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \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 (dimworld == 3)
- * \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;
-
- using std::max;
- const auto dx = b[3] - b[0];
- const auto dy = b[4] - b[1];
- const auto dz = b[5] - b[2];
- const ctype eps = max({dx, dy, dz})*eps_;
- return (b[0] - eps <= point[0] && point[0] <= b[3] + eps &&
- b[1] - eps <= point[1] && point[1] <= b[4] + eps &&
- b[2] - eps <= point[2] && point[2] <= b[5] + eps);
-}
-
-/*!
- * \brief Check whether a point is intersectin a bounding box (dimworld == 2)
- * \param point The point
- * \param b Pointer to bounding box coordinates
- */
-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;
-
- using std::max;
- const auto dx = b[2] - b[0];
- const auto dy = b[3] - b[1];
- const ctype eps = max(dx, dy)*eps_;
- return (b[0] - eps <= point[0] && point[0] <= b[2] + eps &&
- b[1] - eps <= point[1] && point[1] <= b[3] + eps);
-}
-
-/*!
- * \brief Check whether a point is intersectin a bounding box (dimworld == 1)
- * \param point The point
- * \param b Pointer to bounding box coordinates
- */
-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;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Determine if a point intersects an axis-aligned bounding box
- * The bounding box is given by the lower left corner (min) and the upper right corner (max)
- */
-template<class ctype, int dimworld>
-inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point,
- const Dune::FieldVector<ctype, dimworld>& min,
- const Dune::FieldVector<ctype, dimworld>& max)
-{
- std::array<ctype, 2*dimworld> bBox;
- std::copy(min.begin(), min.end(), bBox.begin());
- std::copy(max.begin(), max.end(), bBox.begin()+dimworld);
- return intersectsPointBoundingBox(point, bBox.data());
-}
-
-/*!
- * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 3)
- * \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_*std::max(b[3]-b[0], a[3]-a[0]);
- const ctype eps1 = eps_*std::max(b[4]-b[1], a[4]-a[1]);
- const ctype eps2 = eps_*std::max(b[5]-b[2], a[5]-a[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);
-
-}
-
-/*!
- * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 2)
- * \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 == 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_*std::max(b[2]-b[0], a[2]-a[0]);
- const ctype eps1 = eps_*std::max(b[3]-b[1], a[3]-a[1]);
- return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
- b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
-}
-
-/*!
- * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 1)
- * \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 == 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-7;
- const ctype eps0 = eps_*std::max(b[1]-b[0], a[1]-a[0]);
- return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
-}
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/boundingboxtree.hh>
#endif
diff --git a/dumux/common/geometry/diameter.hh b/dumux/common/geometry/diameter.hh
index 8f4eed920295cb1676755fca486d511e6a3f3aaa..fe4045bec6dd0ead48fcfe6361111dae6247db1b 100644
--- a/dumux/common/geometry/diameter.hh
+++ b/dumux/common/geometry/diameter.hh
@@ -19,31 +19,15 @@
* \ingroup Geometry
* \brief A function to compute a geometry's diameter, i.e.
* the longest distance between points of a geometry
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GEOMETRY_DIAMETER_HH
-#define DUMUX_GEOMETRY_DIAMETER_HH
+#ifndef DUMUX_COMMON_GEOMETRY_DIAMETER_HH
+#define DUMUX_COMMON_GEOMETRY_DIAMETER_HH
-#include <algorithm>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/diameter.hh"
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief Computes the longest distance between points of a geometry
- * \note Useful e.g. to compute the maximum cell diameter of a grid
- */
-template<class Geometry>
-typename Geometry::ctype diameter(const Geometry& geo)
-{
- using std::max;
- typename Geometry::ctype h = 0.0;
- for (std::size_t i = 0; i < geo.corners(); ++i)
- for (std::size_t j = i + 1; j < geo.corners(); ++j)
- h = max(h, (geo.corner(i)-geo.corner(j)).two_norm());
-
- return h;
-}
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/diameter.hh>
#endif
diff --git a/dumux/common/geometry/distance.hh b/dumux/common/geometry/distance.hh
index 47fccf32ea9057aed73b1a2b5ce77c33783d20b5..e87ee09b163d88dc1b41baa9222422bdf23fdc9b 100644
--- a/dumux/common/geometry/distance.hh
+++ b/dumux/common/geometry/distance.hh
@@ -18,164 +18,15 @@
* \file
* \ingroup Geometry
* \brief Helper functions for distance queries
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GEOMETRY_DISTANCE_HH
-#define DUMUX_GEOMETRY_DISTANCE_HH
+#ifndef DUMUX_COMMON_GEOMETRY_DISTANCE_HH
+#define DUMUX_COMMON_GEOMETRY_DISTANCE_HH
-#include <dune/common/fvector.hh>
-#include <dune/geometry/quadraturerules.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/distance.hh"
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief Compute the average distance from a point to a geometry by integration
- */
-template<class Geometry>
-inline typename Geometry::ctype
-averageDistancePointGeometry(const typename Geometry::GlobalCoordinate& p,
- const Geometry& geometry,
- std::size_t integrationOrder = 2)
-{
- typename Geometry::ctype avgDist = 0.0;
- const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
- for (const auto& qp : quad)
- avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position());
- return avgDist/geometry.volume();
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the distance from a point to a line through the points a and b
- */
-template<class Point>
-inline typename Point::value_type
-distancePointLine(const Point& p, const Point& a, const Point& b)
-{
- const auto ab = b - a;
- const auto t = (p - a)*ab/ab.two_norm2();
- auto proj = a;
- proj.axpy(t, ab);
- return (proj - p).two_norm();
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the distance from a point to a line given by a geometry with two corners
- * \note We currently lack the a representation of a line geometry. This convenience function
- * assumes a segment geometry (with two corners) is passed which represents a line geometry.
- */
-template<class Geometry>
-inline typename Geometry::ctype
-distancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
-{
- static_assert(Geometry::mydimension == 1, "Geometry has to be a line");
- const auto& a = geometry.corner(0);
- const auto& b = geometry.corner(1);
- return distancePointLine(p, a, b);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the distance from a point to the segment connecting the points a and b
- */
-template<class Point>
-inline typename Point::value_type
-distancePointSegment(const Point& p, const Point& a, const Point& b)
-{
- const auto ab = b - a;
- auto t = (p - a)*ab;
-
- if (t <= 0.0)
- return (a - p).two_norm();
-
- const auto lengthSq = ab.two_norm2();
- if (t >= lengthSq)
- return (b - p).two_norm();
-
- auto proj = a;
- proj.axpy(t/lengthSq, ab);
- return (proj - p).two_norm();
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the distance from a point to a given segment geometry
- */
-template<class Geometry>
-inline typename Geometry::ctype
-distancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
-{
- static_assert(Geometry::mydimension == 1, "Geometry has to be a segment");
- const auto& a = geometry.corner(0);
- const auto& b = geometry.corner(1);
- return distancePointSegment(p, a, b);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the shortest distance between two points
- */
-template<class ctype, int dimWorld>
-inline ctype distance(const Dune::FieldVector<ctype, dimWorld>& a,
- const Dune::FieldVector<ctype, dimWorld>& b)
-{ return (a-b).two_norm(); }
-
-
-
-namespace Detail {
-
-// helper struct to compute distance between two geometries, specialized below
-template<class Geo1, class Geo2,
- int dimWorld = Geo1::coorddimension,
- int dim1 = Geo1::mydimension, int dim2 = Geo2::mydimension>
-struct GeometryDistance
-{
- static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
- static auto distance(const Geo1& geo1, const Geo2& geo2)
- {
- DUNE_THROW(Dune::NotImplemented, "Geometry distance computation not implemented for dimworld = "
- << dimWorld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
- }
-};
-
-// distance point-point
-template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 0>
-{
- static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
- static auto distance(const Geo1& geo1, const Geo2& geo2)
- { return Dumux::distance(geo1.corner(0), geo2.corner(0)); }
-};
-
-// distance segment-point
-template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 1, 0>
-{
- static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
- static auto distance(const Geo1& geo1, const Geo2& geo2)
- { return distancePointSegment(geo2.corner(0), geo1); }
-};
-
-// distance point-segment
-template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 1>
-{
- static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
- static auto distance(const Geo1& geo1, const Geo2& geo2)
- { return distancePointSegment(geo1.corner(0), geo2); }
-};
-
-} // end namespace Detail
-
-/*!
- * \ingroup Geometry
- * \brief Compute the shortest distance between two geometries
- */
-template<class Geo1, class Geo2>
-inline auto distance(const Geo1& geo1, const Geo2& geo2)
-{ return Detail::GeometryDistance<Geo1, Geo2>::distance(geo1, geo2); }
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/distance.hh>
#endif
diff --git a/dumux/common/geometry/geometricentityset.hh b/dumux/common/geometry/geometricentityset.hh
index 012ace5fccdaefdc3531ec55e49cb379a201a6e1..c17166fe0a99e943365e7e42d8d3a98d51f1ff10 100644
--- a/dumux/common/geometry/geometricentityset.hh
+++ b/dumux/common/geometry/geometricentityset.hh
@@ -18,222 +18,15 @@
* \file
* \ingroup Geometry
* \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
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GEOMETRIC_ENTITY_SET_HH
-#define DUMUX_GEOMETRIC_ENTITY_SET_HH
+#ifndef DUMUX_COMMON_GEOMETRY_GEOMETRIC_ENTITY_SET_HH
+#define DUMUX_COMMON_GEOMETRY_GEOMETRIC_ENTITY_SET_HH
-#include <memory>
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/geometry/multilineargeometry.hh>
-#include <dumux/common/entitymap.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/geometricentityset.hh"
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief An interface for a set of geometric entities based on a GridView
- * \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
- */
-template <class GridView, int codim = 0, class Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>>
-class GridViewGeometricEntitySet
-{
- using EntityMap = Dumux::EntityMap<GridView, codim>;
-public:
- using Entity = typename GridView::template Codim<codim>::Entity;
-
- GridViewGeometricEntitySet(const GridView& gridView)
- : GridViewGeometricEntitySet(gridView, Mapper(gridView, Dune::mcmgLayout(Dune::Codim<codim>())))
- {}
-
- GridViewGeometricEntitySet(const GridView& gridView, const Mapper& mapper)
- : gridView_(gridView)
- , mapper_(mapper)
- , entityMap_(std::make_shared<EntityMap>(gridView.grid(), mapper_))
- {}
-
- GridViewGeometricEntitySet(const GridView& gridView,
- const Mapper& mapper,
- std::shared_ptr<const EntityMap> 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> entityMap_;
-
-};
-
-/*!
- * \ingroup Geometry
- * \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
- */
-template<class GeoType>
-class GeometriesEntitySet
-{
- /*!
- * \brief Wrapper to turn a geometry into a geometric entity
- */
- class EntityWrapper
- {
- public:
- using Geometry = GeoType;
-
- /*!
- * \brief Constructor
- */
- EntityWrapper(const Geometry& geo, const std::size_t index) : geo_(geo), index_(index) {}
-
- /*!
- * \brief Constructor
- */
- EntityWrapper(Geometry&& geo, const std::size_t index) : geo_(std::move(geo)), index_(index) {}
-
- /*!
- * \brief Returns the geometry
- */
- const Geometry& geometry() const
- { return geo_; }
-
- /*!
- * \brief Returns the index of the geometry
- */
- std::size_t index() const
- { return index_; }
-
- private:
- Geometry geo_;
- std::size_t index_;
- };
-
-public:
- using Entity = EntityWrapper;
-
- /*!
- * \brief Constructor for initializer_list
- */
- GeometriesEntitySet(std::initializer_list<typename Entity::Geometry>&& geometries)
- {
- std::size_t index = 0;
- // note: std::initializer_list::begin() returns const T*,
- // thus no moving will be performed and only the copying ctor of
- // EntityWrapper can be called
- for (auto&& g : geometries)
- entities_.emplace_back(g, index++);
- }
-
- /*!
- * \brief Constructor for a vector of geometries
- */
- GeometriesEntitySet(const std::vector<typename Entity::Geometry>& geometries)
- {
- std::size_t index = 0;
- for (auto&& g : geometries)
- entities_.emplace_back(g, index++);
- }
-
- /*!
- * \brief Constructor for a vector of geometries
- */
- GeometriesEntitySet(std::vector<typename Entity::Geometry>&& geometries)
- {
- std::size_t index = 0;
- for (auto&& g : geometries)
- entities_.emplace_back(std::move(g), index++);
- }
-
- /*!
- * \brief The world dimension of the entity set
- */
- enum { dimensionworld = Entity::Geometry::coorddimension };
-
- /*!
- * \brief the coordinate type
- */
- using ctype = typename Entity::Geometry::ctype;
-
- /*!
- * \brief the number of entities in this set
- */
- decltype(auto) size() const
- { return entities_.size(); }
-
- /*!
- * \brief begin iterator to enable range-based for iteration
- */
- decltype(auto) begin() const
- { return entities_.begin(); }
-
- /*!
- * \brief end iterator to enable range-based for iteration
- */
- decltype(auto) end() const
- { return entities_.end(); }
-
- /*!
- * \brief get an entities index
- */
- template<class Entity>
- std::size_t index(const Entity& e) const
- { return e.index(); }
-
- /*!
- * \brief get an entity from an index
- */
- Entity entity(std::size_t index) const
- { return entities_[index]; }
-
-private:
- std::vector<Entity> entities_;
-};
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/geometricentityset.hh>
#endif
diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index 31ce8ea999756f07a1d6c52c94d64c7ca3f377c7..97b350cd78e63dd9b24e3d4a0cb24773226c269b 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -19,1369 +19,15 @@
* \ingroup Geometry
* \brief A class for collision detection of two geometries
* and computation of intersection corners
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GEOMETRY_INTERSECTION_HH
-#define DUMUX_GEOMETRY_INTERSECTION_HH
+#ifndef DUMUX_COMMON_GEOMETRY_INTERSECTION_HH
+#define DUMUX_COMMON_GEOMETRY_INTERSECTION_HH
-#include <tuple>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/geometryintersection.hh"
-#include <dune/common/exceptions.hh>
-#include <dune/common/promotiontraits.hh>
-#include <dune/geometry/multilineargeometry.hh>
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/geometryintersection.hh>
-#include <dumux/common/math.hh>
-#include <dumux/common/geometry/intersectspointgeometry.hh>
-#include <dumux/common/geometry/grahamconvexhull.hh>
-#include <dumux/common/geometry/boundingboxtree.hh>
-
-namespace Dumux {
-namespace IntersectionPolicy {
-
-//! Policy structure for point-like intersections
-template<class ct, int dw>
-struct PointPolicy
-{
- using ctype = ct;
- using Point = Dune::FieldVector<ctype, dw>;
-
- static constexpr int dimWorld = dw;
- static constexpr int dimIntersection = 0;
- using Intersection = Point;
-};
-
-//! Policy structure for segment-like intersections
-template<class ct, int dw>
-struct SegmentPolicy
-{
- using ctype = ct;
- using Point = Dune::FieldVector<ctype, dw>;
- using Segment = std::array<Point, 2>;
-
- static constexpr int dimWorld = dw;
- static constexpr int dimIntersection = 1;
- using Intersection = Segment;
-};
-
-//! Policy structure for polygon-like intersections
-template<class ct, int dw>
-struct PolygonPolicy
-{
- using ctype = ct;
- using Point = Dune::FieldVector<ctype, dw>;
- using Polygon = std::vector<Point>;
-
- static constexpr int dimWorld = dw;
- static constexpr int dimIntersection = 2;
- using Intersection = Polygon;
-};
-
-//! Policy structure for polyhedron-like intersections
-template<class ct, int dw>
-struct PolyhedronPolicy
-{
- using ctype = ct;
- using Point = Dune::FieldVector<ctype, dw>;
- using Polyhedron = std::vector<Point>;
-
- static constexpr int dimWorld = dw;
- static constexpr int dimIntersection = 3;
- using Intersection = Polyhedron;
-};
-
-//! default policy chooser class
-template<class Geometry1, class Geometry2>
-class DefaultPolicyChooser
-{
- static constexpr int dimworld = Geometry1::coorddimension;
- static constexpr int isDim = std::min( int(Geometry1::mydimension), int(Geometry2::mydimension) );
- static_assert(dimworld == int(Geometry2::coorddimension),
- "Geometries must have the same coordinate dimension!");
- static_assert(int(Geometry1::mydimension) <= 3 && int(Geometry2::mydimension) <= 3,
- "Geometries must have dimension 3 or less.");
- using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
-
- using DefaultPolicies = std::tuple<PointPolicy<ctype, dimworld>,
- SegmentPolicy<ctype, dimworld>,
- PolygonPolicy<ctype, dimworld>,
- PolyhedronPolicy<ctype, dimworld>>;
-public:
- using type = std::tuple_element_t<isDim, DefaultPolicies>;
-};
-
-//! Helper alias to define the default intersection policy
-template<class Geometry1, class Geometry2>
-using DefaultPolicy = typename DefaultPolicyChooser<Geometry1, Geometry2>::type;
-
-} // end namespace IntersectionPolicy
-
-namespace Detail {
-
-/*!
- * \ingroup Geometry
- * \brief Algorithm to find segment-like intersections of a polgon/polyhedron with a
- * segment. The result is stored in the form of the local coordinates tfirst
- * and tlast on the segment geo1.
- * \param geo1 the first geometry
- * \param geo2 the second geometry (dimGeo2 < dimGeo1)
- * \param baseEps the base epsilon used for floating point comparisons
- * \param tfirst stores the local coordinate of beginning of intersection segment
- * \param tlast stores the local coordinate of the end of intersection segment
- * \param getFacetCornerIndices Function to obtain the facet corner indices on geo1
- * \param computeNormal Function to obtain the normal vector on a facet on geo1
- */
-template< class Geo1, class Geo2, class ctype,
- class GetFacetCornerIndices, class ComputeNormalFunction >
-bool computeSegmentIntersection(const Geo1& geo1, const Geo2& geo2, ctype baseEps,
- ctype& tfirst, ctype& tlast,
- const GetFacetCornerIndices& getFacetCornerIndices,
- const ComputeNormalFunction& computeNormal)
-{
- static_assert(int(Geo2::mydimension) == 1, "Geometry2 must be a segment");
- static_assert(int(Geo1::mydimension) > int(Geo2::mydimension),
- "Dimension of Geometry1 must be higher than that of Geometry2");
-
- const auto a = geo2.corner(0);
- const auto b = geo2.corner(1);
- const auto d = b - a;
-
- // The initial interval is the whole segment.
- // Afterwards we start clipping the interval
- // at the intersections with the facets of geo1
- tfirst = 0.0;
- tlast = 1.0;
-
- const auto facets = getFacetCornerIndices(geo1);
- for (const auto& f : facets)
- {
- const auto n = computeNormal(f);
-
- const auto c0 = geo1.corner(f[0]);
- const ctype denom = n*d;
- const ctype dist = n*(a-c0);
-
- // use first edge of the facet to scale eps
- const auto edge1 = geo1.corner(f[1]) - geo1.corner(f[0]);
- const ctype eps = baseEps*edge1.two_norm();
-
- // if denominator is zero the segment in parallel to the edge.
- // If the distance is positive there is no intersection
- using std::abs;
- if (abs(denom) < eps)
- {
- if (dist > eps)
- return false;
- }
- else // not parallel: compute line-line intersection
- {
- const ctype t = -dist / denom;
- // if entering half space cut tfirst if t is larger
- using std::signbit;
- if (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 - baseEps)
- return false;
- }
- }
-
- // If we made it until here an intersection exists.
- return true;
-}
-
-} // end namespace Detail
-
-/*!
- * \ingroup Geometry
- * \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 Geometry1, class Geometry2,
- class Policy = IntersectionPolicy::DefaultPolicy<Geometry1, Geometry2>,
- int dimworld = Geometry1::coorddimension,
- int dim1 = Geometry1::mydimension,
- int dim2 = Geometry2::mydimension>
-class GeometryIntersection
-{
- static constexpr int dimWorld = Policy::dimWorld;
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
- //! Determine if the two geometries intersect and compute the intersection geometry
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- 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);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for segment--segment intersection in 2d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 1>
-{
- enum { dimworld = 2 };
- enum { dim1 = 1 };
- enum { dim2 = 1 };
-
- // base epsilon for floating point comparisons
- static constexpr typename Policy::ctype eps_ = 1.5e-7;
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
- /*!
- * \brief Colliding two segments
- * \param geo1/geo2 The geometries to intersect
- * \param intersection The intersection point
- * \note This overload is used when point-like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- const auto v1 = geo1.corner(1) - geo1.corner(0);
- const auto v2 = geo2.corner(1) - geo2.corner(0);
- const auto ac = geo2.corner(0) - geo1.corner(0);
-
- auto n2 = Point({-1.0*v2[1], v2[0]});
- n2 /= n2.two_norm();
-
- // compute distance of first corner on geo2 to line1
- const auto dist12 = n2*ac;
-
- // first check parallel segments
- using std::abs;
- using std::sqrt;
-
- const auto v1Norm2 = v1.two_norm2();
- const auto eps = eps_*sqrt(v1Norm2);
- const auto eps2 = eps_*v1Norm2;
-
- const auto sp = n2*v1;
- if (abs(sp) < eps)
- {
- if (abs(dist12) > eps)
- return false;
-
- // point intersection can only be one of the corners
- if ( (v1*v2) < 0.0 )
- {
- if ( ac.two_norm2() < eps2 )
- { intersection = geo2.corner(0); return true; }
-
- if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 )
- { intersection = geo2.corner(1); return true; }
- }
- else
- {
- if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 )
- { intersection = geo2.corner(1); return true; }
-
- if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 )
- { intersection = geo2.corner(0); return true; }
- }
-
- // no intersection
- return false;
- }
-
- // intersection point on v1 in local coords
- const auto t1 = dist12 / sp;
-
- // check if the local coords are valid
- if (t1 < -1.0*eps_ || t1 > 1.0 + eps_)
- return false;
-
- // compute global coordinates
- auto isPoint = geo1.global(t1);
-
- // check if point is in bounding box of v2
- const auto c = geo2.corner(0);
- const auto d = geo2.corner(1);
-
- using std::min; using std::max;
- std::array<ctype, 4> bBox({ min(c[0], d[0]), min(c[1], d[1]), max(c[0], d[0]), max(c[1], d[1]) });
-
- if ( intersectsPointBoundingBox(isPoint, bBox.data()) )
- {
- intersection = std::move(isPoint);
- return true;
- }
-
- return false;
- }
-
- /*!
- * \brief Colliding two segments
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store corners of intersection segment
- * \note this overload is used when segment-like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- const auto& a = geo1.corner(0);
- const auto& b = geo1.corner(1);
- const auto ab = b - a;
-
- const auto& p = geo2.corner(0);
- const auto& q = geo2.corner(1);
- const auto pq = q - p;
- Dune::FieldVector<ctype, dimworld> n2{-pq[1], pq[0]};
-
- using std::max;
- const auto abNorm2 = ab.two_norm2();
- const auto pqNorm2 = pq.two_norm2();
- const auto eps2 = eps_*max(abNorm2, pqNorm2);
-
- // non-parallel segments do not intersect in a segment.
- using std::abs;
- if (abs(n2*ab) > eps2)
- return false;
-
- // check if the segments lie on the same line
- const auto ap = p - a;
- if (abs(ap*n2) > eps2)
- return false;
-
- // compute scaled local coordinates of corner 1/2 of segment2 on segment1
- auto t1 = ab*ap;
- auto t2 = ab*(q - a);
-
- using std::swap;
- if (t1 > t2)
- swap(t1, t2);
-
- using std::clamp;
- t1 = clamp(t1, 0.0, abNorm2);
- t2 = clamp(t2, 0.0, abNorm2);
-
- if (abs(t2-t1) < eps2)
- return false;
-
- intersection = Intersection({geo1.global(t1/abNorm2), geo1.global(t2/abNorm2)});
- return true;
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polygon--segment intersection in 2d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 1>
-{
- enum { dimworld = 2 };
- enum { dim1 = 2 };
- enum { dim2 = 1 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding segment and convex polygon
- * \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.
- * \note This overload is used when segment-like intersections are seeked.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- ctype tfirst, tlast;
- if (intersect_(geo1, geo2, tfirst, tlast))
- {
- // the intersection exists. Export the intersections geometry now:
- // s(t) = a + t(b-a) in [tfirst, tlast]
- intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
- return true;
- }
-
- return false;
- }
-
- /*!
- * \brief Colliding segment and convex polygon
- * \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.
- * \note This overload is used when point-like intersections are seeked.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename P::Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- ctype tfirst, tlast;
- if (!intersect_(geo1, geo2, tfirst, tlast))
- {
- // if start & end point are same, the intersection is a point
- using std::abs;
- if ( tfirst > tlast - eps_ )
- {
- intersection = Intersection(geo2.global(tfirst));
- return true;
- }
- }
-
- return false;
- }
-
-private:
- /*!
- * \brief Obtain local coordinates of start/end point of the intersecting segment.
- */
- static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
- {
- // lambda to obtain the facet corners on geo1
- auto getFacetCorners = [] (const Geometry1& geo1)
- {
- std::vector< std::array<int, 2> > facetCorners;
- switch (geo1.corners())
- {
- case 4: // quadrilateral
- facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
- break;
- case 3: // triangle
- facetCorners = {{1, 0}, {0, 2}, {2, 1}};
- break;
- default:
- DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
- << geo1.type() << " with "<< geo1.corners() << " corners.");
- }
-
- return facetCorners;
- };
-
- // lambda to obtain the normal on a facet
- const auto center1 = geo1.center();
- auto computeNormal = [&geo1, ¢er1] (const std::array<int, 2>& facetCorners)
- {
- const auto c0 = geo1.corner(facetCorners[0]);
- const auto c1 = geo1.corner(facetCorners[1]);
- const auto edge = c1 - c0;
-
- Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]});
- n /= n.two_norm();
-
- // orientation of the element is not known
- // make sure normal points outwards of element
- if ( n*(center1-c0) > 0.0 )
- n *= -1.0;
-
- return n;
- };
-
- return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for segment--polygon intersection in 2d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 2>
-: public GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1>
-{
- using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 2, 2, 1>;
-
-public:
- /*!
- * \brief Colliding segment and convex polygon
- * \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.
- * \note This forwards to the polygon-segment specialization with swapped arguments.
- */
- template<class P = Policy>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
- {
- return Base::intersection(geo2, geo1, intersection);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polygon--polygon intersection in 2d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 2>
-{
- enum { dimworld = 2 };
- enum { dim1 = 2 };
- enum { dim2 = 2 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding two polygons
- * \note First we find the vertex candidates for the intersection region as follows:
- * Add polygon vertices that are inside the other polygon
- * Add intersections of polygon edges
- * Remove duplicate points from the list
- * Compute the convex hull polygon
- * Return a triangulation of that polygon as intersection
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store the corner points of the polygon (as convex hull)
- * \note This overload is used when polygon like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- // the candidate intersection points
- std::vector<Point> points; points.reserve(6);
-
- // add polygon1 corners that are inside polygon2
- for (int i = 0; i < geo1.corners(); ++i)
- if (intersectsPointGeometry(geo1.corner(i), geo2))
- points.emplace_back(geo1.corner(i));
-
- const auto numPoints1 = points.size();
- if (numPoints1 != geo1.corners())
- {
- // add polygon2 corners that are inside polygon1
- for (int i = 0; i < geo2.corners(); ++i)
- if (intersectsPointGeometry(geo2.corner(i), geo1))
- points.emplace_back(geo2.corner(i));
-
- if (points.empty())
- return false;
-
- if (points.size() - numPoints1 != geo2.corners())
- {
- const auto refElement1 = referenceElement(geo1);
- const auto refElement2 = referenceElement(geo2);
-
- // add intersections of edges
- using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
- using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
- for (int i = 0; i < refElement1.size(dim1-1); ++i)
- {
- const auto localEdgeGeom1 = refElement1.template geometry<dim1-1>(i);
- const auto edge1 = SegGeometry( Dune::GeometryTypes::line,
- std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)),
- geo1.global(localEdgeGeom1.corner(1))} ));
-
- for (int j = 0; j < refElement2.size(dim2-1); ++j)
- {
- const auto localEdgeGeom2 = refElement2.template geometry<dim2-1>(j);
- const auto edge2 = SegGeometry( Dune::GeometryTypes::line,
- std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)),
- geo2.global(localEdgeGeom2.corner(1))} ));
-
- using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>;
- typename EdgeTest::Intersection edgeIntersection;
- if (EdgeTest::intersection(edge1, edge2, edgeIntersection))
- points.emplace_back(edgeIntersection);
- }
- }
- }
- }
-
- if (points.empty())
- return false;
-
- // remove duplicates
- const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
- std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
- {
- using std::abs;
- return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : a[1] < b[1]);
- });
-
- auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
- {
- return (b-a).two_norm() < eps;
- });
-
- points.erase(removeIt, points.end());
-
- // return false if we don't have at least three unique points
- if (points.size() < 3)
- return false;
-
- // intersection polygon is convex hull of above points
- intersection = grahamConvexHull<2>(points);
- assert(!intersection.empty());
- return true;
- }
-
- /*!
- * \brief Colliding two polygons
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store the corners of intersection segment
- * \note this overload is used when segment-like intersections are seeked
- * \todo implement this query
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for segment-like intersections");
- }
-
- /*!
- * \brief Colliding two polygons
- * \param geo1/geo2 The geometries to intersect
- * \param intersection The intersection point
- * \note this overload is used when point-like intersections are seeked
- * \todo implement this query
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for touching points");
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polyhedron--segment intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 1>
-{
- enum { dimworld = 3 };
- enum { dim1 = 3 };
- enum { dim2 = 1 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding segment and convex polyhedron
- * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
- * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
- * 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.
- * \note This overload is used when segment-like intersections are seeked.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- ctype tfirst, tlast;
- if (intersect_(geo1, geo2, tfirst, tlast))
- {
- // Intersection exists. Export the intersections geometry now:
- // s(t) = a + t(b-a) in [tfirst, tlast]
- intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
- return true;
- }
-
- return false;
- }
-
- /*!
- * \brief Colliding segment and convex polyhedron
- * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
- * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
- * 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.
- * \note This overload is used when point-like intersections are seeked.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- ctype tfirst, tlast;
- if (intersect_(geo1, geo2, tfirst, tlast))
- {
- // if start & end point are same, the intersection is a point
- using std::abs;
- if ( abs(tfirst - tlast) < eps_ )
- {
- intersection = Intersection(geo2.global(tfirst));
- return true;
- }
- }
-
- return false;
- }
-
-private:
- /*!
- * \brief Obtain local coordinates of start/end point of the intersecting segment.
- */
- static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
- {
- // lambda to obtain facet corners on geo1
- auto getFacetCorners = [] (const Geometry1& geo1)
- {
- std::vector< std::vector<int> > facetCorners;
- // 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
- facetCorners = {{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
- facetCorners = {{1, 0, 2}, {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 facetCorners;
- };
-
- // lambda to obtain the normal on a facet
- auto computeNormal = [&geo1] (const std::vector<int>& facetCorners)
- {
- const auto v0 = geo1.corner(facetCorners[1]) - geo1.corner(facetCorners[0]);
- const auto v1 = geo1.corner(facetCorners[2]) - geo1.corner(facetCorners[0]);
-
- auto n = crossProduct(v0, v1);
- n /= n.two_norm();
-
- return n;
- };
-
- return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for segment--polyhedron intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 3>
-: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1>
-{
- using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 1>;
-public:
- /*!
- * \brief Colliding segment and convex polyhedron
- * \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.
- * \note This forwards to the polyhedron-segment specialization with swapped arguments.
- */
- template<class P = Policy>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
- {
- return Base::intersection(geo2, geo1, intersection);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polyhedron--polygon intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 2>
-{
- enum { dimworld = 3 };
- enum { dim1 = 3 };
- enum { dim2 = 2 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding polygon and convex polyhedron
- * \note First we find the vertex candidates for the intersection region as follows:
- * Add triangle vertices that are inside the tetrahedron
- * Add tetrahedron vertices that are inside the triangle
- * Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests)
- * Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests)
- * Remove duplicate points from the list
- * Compute the convex hull polygon
- * Return a triangulation of that polygon as intersection
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store the corner points of the polygon (as convex hull)
- * \note This overload is used when polygon like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- // the candidate intersection points
- std::vector<Point> points; points.reserve(10);
-
- // add 3d geometry corners that are inside the 2d geometry
- for (int i = 0; i < geo1.corners(); ++i)
- if (intersectsPointGeometry(geo1.corner(i), geo2))
- points.emplace_back(geo1.corner(i));
-
- // add 2d geometry corners that are inside the 3d geometry
- for (int i = 0; i < geo2.corners(); ++i)
- if (intersectsPointGeometry(geo2.corner(i), geo1))
- points.emplace_back(geo2.corner(i));
-
- // get some geometry types
- using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>;
- using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
-
- const auto refElement1 = referenceElement(geo1);
- const auto refElement2 = referenceElement(geo2);
-
- // intersection policy for point-like intersections (used later)
- using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
-
- // add intersection points of all polyhedron edges (codim dim-1) with the polygon
- for (int i = 0; i < refElement1.size(dim1-1); ++i)
- {
- const auto localEdgeGeom = refElement1.template geometry<dim1-1>(i);
- const auto p = geo1.global(localEdgeGeom.corner(0));
- const auto q = geo1.global(localEdgeGeom.corner(1));
- const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
-
- using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>;
- typename PolySegTest::Intersection polySegIntersection;
- if (PolySegTest::intersection(geo2, segGeo, polySegIntersection))
- points.emplace_back(polySegIntersection);
- }
-
- // add intersection points of all polygon faces (codim 1) with the polyhedron faces
- for (int i = 0; i < refElement1.size(1); ++i)
- {
- const auto faceGeo = [&]()
- {
- const auto localFaceGeo = refElement1.template geometry<1>(i);
- if (localFaceGeo.corners() == 4)
- {
- const auto a = geo1.global(localFaceGeo.corner(0));
- const auto b = geo1.global(localFaceGeo.corner(1));
- const auto c = geo1.global(localFaceGeo.corner(2));
- const auto d = geo1.global(localFaceGeo.corner(3));
-
- return PolyhedronFaceGeometry(Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d});
- }
- else
- {
- const auto a = geo1.global(localFaceGeo.corner(0));
- const auto b = geo1.global(localFaceGeo.corner(1));
- const auto c = geo1.global(localFaceGeo.corner(2));
-
- return PolyhedronFaceGeometry(Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c});
- }
- }();
-
- for (int j = 0; j < refElement2.size(1); ++j)
- {
- const auto localEdgeGeom = refElement2.template geometry<1>(j);
- const auto p = geo2.global(localEdgeGeom.corner(0));
- const auto q = geo2.global(localEdgeGeom.corner(1));
-
- const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
-
- using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>;
- typename PolySegTest::Intersection polySegIntersection;
- if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection))
- points.emplace_back(polySegIntersection);
- }
- }
-
- // return if no intersection points were found
- if (points.empty()) return false;
-
- // remove duplicates
- const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
- std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
- {
- using std::abs;
- return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
- });
-
- auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
- {
- return (b-a).two_norm() < eps;
- });
-
- points.erase(removeIt, points.end());
-
- // return false if we don't have more than three unique points
- if (points.size() < 3) return false;
-
- // intersection polygon is convex hull of above points
- intersection = grahamConvexHull<2>(points);
- assert(!intersection.empty());
-
- return true;
- }
-
- /*!
- * \brief Colliding segment and convex polyhedron
- * \note First we find the vertex candidates for the intersection region as follows:
- * Add triangle vertices that are inside the tetrahedron
- * Add tetrahedron vertices that are inside the triangle
- * Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests)
- * Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests)
- * Remove duplicate points from the list
- * Compute the convex hull polygon
- * Return a triangulation of that polygon as intersection
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store the intersection result
- * \todo implement overloads for segment or point-like intersections
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection != 2, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- DUNE_THROW(Dune::NotImplemented, "Polyhedron-polygon intersection detection only implemented for polygon-like intersections");
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polygon--polyhedron intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 3>
-: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2>
-{
- using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 3, 2>;
-public:
- /*!
- * \brief Colliding polygon and convex polyhedron
- * \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.
- * \note This forwards to the polyhedron-polygon specialization with swapped arguments.
- */
- template<class P = Policy>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
- {
- return Base::intersection(geo2, geo1, intersection);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for polygon--segment intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 1>
-{
- enum { dimworld = 3 };
- enum { dim1 = 2 };
- enum { dim2 = 1 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding segment and convex polyhedron
- * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
- * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.3.6)
- * \param geo1/geo2 The geometries to intersect
- * \param intersection If the geometries collide, is holds the corner points of
- * the intersection object in global coordinates.
- * \note This overload is used when point-like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- ctype tfirst, tlast;
- if (intersect_(geo1, geo2, tfirst, tlast))
- {
- // the intersection exists. Export the intersections geometry now:
- // s(t) = a + t(b-a) in [tfirst, tlast]
- intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
- return true;
- }
-
- return false;
- }
-
- /*!
- * \brief Colliding segment and convex polyhedron
- * \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
- * published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.3.6)
- * \param geo1/geo2 The geometries to intersect
- * \param is If the geometries collide, is holds the corner points of
- * the intersection object in global coordinates.
- * \note This overload is used when point-like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& is)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- const auto p = geo2.corner(0);
- const auto q = geo2.corner(1);
-
- const auto a = geo1.corner(0);
- const auto b = geo1.corner(1);
- const auto c = geo1.corner(2);
-
- if (geo1.corners() == 3)
- return intersect_<Policy>(a, b, c, p, q, is);
-
- else if (geo1.corners() == 4)
- {
- Intersection is1, is2;
- bool hasSegment1, hasSegment2;
-
- const auto d = geo1.corner(3);
- const bool intersects1 = intersect_<Policy>(a, b, d, p, q, is1, hasSegment1);
- const bool intersects2 = intersect_<Policy>(a, d, c, p, q, is2, hasSegment2);
-
- if (hasSegment1 || hasSegment2)
- return false;
-
- if (intersects1) { is = is1; return true; }
- if (intersects2) { is = is2; return true; }
-
- return false;
- }
-
- else
- DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
- << geo1.type() << ", "<< geo1.corners() << " corners.");
- }
-
-private:
- /*!
- * \brief Obtain local coordinates of start/end point of the intersecting segment.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
- {
- // lambda to obtain the facet corners on geo1
- auto getFacetCorners = [] (const Geometry1& geo1)
- {
- std::vector< std::array<int, 2> > facetCorners;
- switch (geo1.corners())
- {
- case 4: // quadrilateral
- facetCorners = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
- break;
- case 3: // triangle
- facetCorners = {{1, 0}, {0, 2}, {2, 1}};
- break;
- default:
- DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
- << geo1.type() << " with "<< geo1.corners() << " corners.");
- }
-
- return facetCorners;
- };
-
- const auto center1 = geo1.center();
- const auto normal1 = crossProduct(geo1.corner(1) - geo1.corner(0),
- geo1.corner(2) - geo1.corner(0));
-
- // lambda to obtain the normal on a facet
- auto computeNormal = [¢er1, &normal1, &geo1] (const std::array<int, 2>& facetCorners)
- {
- const auto c0 = geo1.corner(facetCorners[0]);
- const auto c1 = geo1.corner(facetCorners[1]);
- const auto edge = c1 - c0;
-
- Dune::FieldVector<ctype, dimworld> n = crossProduct(edge, normal1);
- n /= n.two_norm();
-
- // orientation of the element is not known
- // make sure normal points outwards of element
- if ( n*(center1-c0) > 0.0 )
- n *= -1.0;
-
- return n;
- };
-
- return Detail::computeSegmentIntersection(geo1, geo2, eps_, tfirst, tlast, getFacetCorners, computeNormal);
- }
-
- /*!
- * \brief triangle--segment point-like intersection with points as input.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersect_(const Point& a, const Point& b, const Point& c,
- const Point& p, const Point& q,
- Intersection& is)
- {
- bool hasSegment;
- return intersect_(a, b, c, p, q, is, hasSegment);
- }
-
- /*!
- * \brief triangle--segment point-like intersection with points as input. Also
- * stores if a segment-like intersection was found in the provided boolean.
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersect_(const Point& a, const Point& b, const Point& c,
- const Point& p, const Point& q,
- Intersection& is, bool& hasSegmentIs)
- {
- hasSegmentIs = false;
-
- const auto ab = b - a;
- const auto ac = c - a;
- const auto qp = p - q;
-
- // compute the triangle normal that defines the triangle plane
- const auto normal = crossProduct(ab, ac);
-
- // compute the denominator
- // if denom is 0 the segment is parallel and we can return
- const auto denom = normal*qp;
- const auto abnorm2 = ab.two_norm2();
- const auto eps = eps_*abnorm2*qp.two_norm();
- using std::abs;
- if (abs(denom) < eps)
- {
- const auto pa = a - p;
- const auto denom2 = normal*pa;
- if (abs(denom2) > eps_*pa.two_norm()*abnorm2)
- return false;
-
- // if we get here, we are in-plane. Check if a
- // segment-like intersection with geometry 1 exists.
- using SegmentPolicy = typename IntersectionPolicy::SegmentPolicy<ctype, dimworld>;
- using Triangle = Dune::AffineGeometry<ctype, 2, dimworld>;
- using Segment = Dune::AffineGeometry<ctype, 1, dimworld>;
- using SegmentIntersectionAlgorithm = GeometryIntersection<Triangle, Segment, SegmentPolicy>;
- using SegmentIntersectionType = typename SegmentIntersectionAlgorithm::Intersection;
- SegmentIntersectionType segmentIs;
-
- Triangle triangle(Dune::GeometryTypes::simplex(2), std::array<Point, 3>({a, b, c}));
- Segment segment(Dune::GeometryTypes::simplex(1), std::array<Point, 2>({p, q}));
- if (SegmentIntersectionAlgorithm::intersection(triangle, segment, segmentIs))
- {
- hasSegmentIs = true;
- return false;
- }
-
- // We are in-plane. A point-like
- // intersection can only be on the edges
- for (const auto& ip : {p, q})
- {
- if ( intersectsPointSimplex(ip, a, b)
- || intersectsPointSimplex(ip, b, c)
- || intersectsPointSimplex(ip, c, a) )
- {
- is = ip;
- return true;
- }
- }
-
- return false;
- }
-
- // compute intersection t value of pq with plane of triangle.
- // a segment intersects if and only if 0 <= t <= 1.
- const auto ap = p - a;
- const auto t = (ap*normal)/denom;
- if (t < 0.0 - eps_) return false;
- if (t > 1.0 + eps_) return false;
-
- // compute the barycentric coordinates and check if the intersection point
- // is inside the bounds of the triangle
- const auto e = crossProduct(qp, ap);
- const auto v = (ac*e)/denom;
- if (v < -eps_ || v > 1.0 + eps_) return false;
- const auto w = -(ab*e)/denom;
- if (w < -eps_ || v + w > 1.0 + eps_) return false;
-
- // Now we are sure there is an intersection points
- // Perform delayed division compute the last barycentric coordinate component
- const auto u = 1.0 - v - w;
-
- Point ip(0.0);
- ip.axpy(u, a);
- ip.axpy(v, b);
- ip.axpy(w, c);
- is = ip;
- return true;
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for segment--polygon intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 2>
-: public GeometryIntersection<Geometry2, Geometry1, Policy, 3, 2, 1>
-{
- using Base = GeometryIntersection<Geometry2, Geometry1, Policy, 3, 2, 1>;
-public:
- /*!
- * \brief Colliding segment and convex polygon
- * \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.
- * \note This forwards to the polyhedron-polygon specialization with swapped arguments.
- */
- template<class P = Policy>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename Base::Intersection& intersection)
- {
- return Base::intersection(geo2, geo1, intersection);
- }
-};
-
-/*!
- * \ingroup Geometry
- * \brief A class for segment--segment intersection in 3d space
- */
-template <class Geometry1, class Geometry2, class Policy>
-class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 1, 1>
-{
- enum { dimworld = 3 };
- enum { dim1 = 1 };
- enum { dim2 = 1 };
-
-public:
- using ctype = typename Policy::ctype;
- using Point = typename Policy::Point;
- using Intersection = typename Policy::Intersection;
-
-private:
- static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-
-public:
- /*!
- * \brief Colliding two segments
- * \param geo1/geo2 The geometries to intersect
- * \param intersection Container to store corners of intersection segment
- * \note this overload is used when point-like intersections are seeked
- * \todo implement this query
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- DUNE_THROW(Dune::NotImplemented, "segment-segment intersection detection for point-like intersections");
- }
-
- /*!
- * \brief Colliding two segments in 3D
- * \param geo1/geo2 The geometries to intersect
- * \param intersection If the geometries collide, is holds the corner points of
- * the intersection object in global coordinates.
- * \note This overload is used when segment-like intersections are seeked
- */
- template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
- {
- static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
-
- const auto& a = geo1.corner(0);
- const auto& b = geo1.corner(1);
- const auto ab = b-a;
-
- const auto& p = geo2.corner(0);
- const auto& q = geo2.corner(1);
- const auto pq = q-p;
-
- const auto abNorm2 = ab.two_norm2();
- const auto pqNorm2 = pq.two_norm2();
-
- using std::max;
- const auto eps2 = eps_*max(abNorm2, pqNorm2);
-
- // if the segment are not parallel there is no segment intersection
- using std::abs;
- if (crossProduct(ab, pq).two_norm2() > eps2*eps2)
- return false;
-
- const auto ap = (p-a);
- const auto aq = (q-a);
-
- // points have to be colinear
- if (crossProduct(ap, aq).two_norm2() > eps2*eps2)
- return false;
-
- // scaled local coordinates
- // we save the division for the normalization for now
- // and do it in the very end if we find an intersection
- auto tp = ap*ab;
- auto tq = aq*ab;
-
- // make sure they are sorted
- using std::swap;
- if (tp > tq)
- swap(tp, tq);
-
- using std::clamp;
- tp = clamp(tp, 0.0, abNorm2);
- tq = clamp(tq, 0.0, abNorm2);
-
- if (abs(tp-tq) < eps2)
- return false;
-
- intersection = Intersection({geo1.global(tp/abNorm2), geo1.global(tq/abNorm2)});
- return true;
- }
-};
-
-} // end namespace Dumux
-
-# endif
+#endif
diff --git a/dumux/common/geometry/grahamconvexhull.hh b/dumux/common/geometry/grahamconvexhull.hh
index 3ed273d9e8a98772345d5d41ea3b7d5783f89830..b1a3155bfe7b4b89ce44fe2e0d18277b61f5f830 100644
--- a/dumux/common/geometry/grahamconvexhull.hh
+++ b/dumux/common/geometry/grahamconvexhull.hh
@@ -19,200 +19,15 @@
* \ingroup Geometry
* \brief A function to compute the convex hull of a point cloud
* and a function to triangulate the polygon area spanned by the convex hull
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GRAHAM_CONVEX_HULL_HH
-#define DUMUX_GRAHAM_CONVEX_HULL_HH
+#ifndef DUMUX_COMMON_GEOMETRY_GRAHAM_CONVEX_HULL_HH
+#define DUMUX_COMMON_GEOMETRY_GRAHAM_CONVEX_HULL_HH
-#include <vector>
-#include <array>
-#include <algorithm>
-#include <iterator>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/grahamconvexhull.hh"
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/grahamconvexhull.hh>
-#include <dumux/common/math.hh>
-#include <dumux/common/geometry/triangulation.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief Returns the orientation of a sequence a-->b-->c in one plane (defined by normal vector)
- * \return -1 if a-->b-->c forms a counter-clockwise turn (given the normal vector)
- * +1 for a clockwise turn,
- * 0 if they are on one line (colinear)
- */
-template<class ctype>
-int getOrientation(const Dune::FieldVector<ctype, 3>& a,
- const Dune::FieldVector<ctype, 3>& b,
- const Dune::FieldVector<ctype, 3>& c,
- const Dune::FieldVector<ctype, 3>& normal)
-{
- const auto d = b-a;
- const auto e = c-b;
- const auto f = Dumux::crossProduct(d, e);
- const auto area = f*normal;
- return Dumux::sign(-area);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the points making up the convex hull around the given set of unordered points
- * \note We assume that all points are coplanar and there are no indentical points in the list
- * \note This algorithm changes the order of the given points a bit
- * as they are unordered anyway this shouldn't matter too much
- */
-template<int dim, class ctype,
- std::enable_if_t<(dim==2), int> = 0>
-std::vector<Dune::FieldVector<ctype, 3>>
-grahamConvexHullImpl(std::vector<Dune::FieldVector<ctype, 3>>& points)
-{
- using Point = Dune::FieldVector<ctype, 3>;
- std::vector<Point> convexHull;
-
- // return empty convex hull
- if (points.size() < 3)
- return convexHull;
-
- // return the points (already just one triangle)
- if (points.size() == 3)
- return points;
-
- // try to compute the normal vector of the plane
- const auto a = points[1] - points[0];
- auto b = points[2] - points[0];
- auto normal = Dumux::crossProduct(a, b);
-
- // make sure the normal vector has non-zero length
- std::size_t k = 2;
- auto norm = normal.two_norm();
- while (norm == 0.0 && k < points.size()-1)
- {
- b = points[++k];
- normal = Dumux::crossProduct(a, b);
- norm = normal.two_norm();
- }
-
- // if all given points are colinear -> return empty convex hull
- if (norm == 0.0)
- return convexHull;
-
- using std::sqrt;
- const auto eps = 1e-7*sqrt(norm);
- normal /= norm;
-
- // find the element with the smallest x coordinate (if x is the same, smallest y coordinate, and so on...)
- auto minIt = std::min_element(points.begin(), points.end(), [&eps](const auto& a, const auto& b)
- {
- using std::abs;
- return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
- });
-
- // swap the smallest element to the front
- std::iter_swap(minIt, points.begin());
-
- // choose the first (min element) as the pivot point
- // sort in counter-clockwise order around pivot point
- const auto pivot = points[0];
- std::sort(points.begin()+1, points.end(), [&](const auto& a, const auto& b)
- {
- const int order = getOrientation(pivot, a, b, normal);
- if (order == 0)
- return (a-pivot).two_norm() < (b-pivot).two_norm();
- else
- return (order == -1);
- });
-
- // push the first three points
- convexHull.reserve(50);
- convexHull.push_back(points[0]);
- convexHull.push_back(points[1]);
- convexHull.push_back(points[2]);
-
- // This is the heart of the algorithm
- // pop_back until the last point in the queue forms a counter-clockwise oriented line
- // with the next two vertices. Then add these points to the queue.
- for (std::size_t i = 3; i < points.size(); ++i)
- {
- Point p = convexHull.back();
- convexHull.pop_back();
- // keep popping until the orientation a->b->currentp is counter-clockwise
- while (getOrientation(convexHull.back(), p, points[i], normal) != -1)
- {
- // make sure the queue doesn't get empty
- if (convexHull.size() == 1)
- {
- // before we reach i=size-1 there has to be a good candidate
- // as not all points are colinear (a non-zero plane normal exists)
- assert(i < points.size()-1);
- p = points[i++];
- }
- else
- {
- p = convexHull.back();
- convexHull.pop_back();
- }
- }
-
- // add back the last popped point and this point
- convexHull.emplace_back(std::move(p));
- convexHull.push_back(points[i]);
- }
-
- return convexHull;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the points making up the convex hull around the given set of unordered points
- * \note This is the specialization for 2d space. Here, we make use of the generic implementation
- * for the case of coplanar points in 3d space (a more efficient implementation could be provided).
- */
-template<int dim, class ctype,
- std::enable_if_t<(dim==2), int> = 0>
-std::vector<Dune::FieldVector<ctype, 2>>
-grahamConvexHullImpl(const std::vector<Dune::FieldVector<ctype, 2>>& points)
-{
- std::vector<Dune::FieldVector<ctype, 3>> points3D;
- points3D.reserve(points.size());
- std::transform(points.begin(), points.end(), std::back_inserter(points3D),
- [](const auto& p) { return Dune::FieldVector<ctype, 3>({p[0], p[1], 0.0}); });
-
- const auto result3D = grahamConvexHullImpl<2>(points3D);
-
- std::vector<Dune::FieldVector<ctype, 2>> result2D;
- result2D.reserve(result3D.size());
- std::transform(result3D.begin(), result3D.end(), std::back_inserter(result2D),
- [](const auto& p) { return Dune::FieldVector<ctype, 2>({p[0], p[1]}); });
-
- return result2D;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the points making up the convex hull around the given set of unordered points
- * \note We assume that all points are coplanar and there are no indentical points in the list
- */
-template<int dim, class ctype, int dimWorld>
-std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
-{
- return grahamConvexHullImpl<dim>(points);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the points making up the convex hull around the given set of unordered points
- * \note We assume that all points are coplanar and there are no indentical points in the list
- * \note This is the overload if we are not allowed to write into the given points vector
- */
-template<int dim, class ctype, int dimWorld>
-std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(const std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
-{
- auto copyPoints = points;
- return grahamConvexHullImpl<dim>(copyPoints);
-}
-
-} // end namespace Dumux
-
-# endif
+#endif
diff --git a/dumux/common/geometry/intersectingentities.hh b/dumux/common/geometry/intersectingentities.hh
index 15378819faf88a687851de7877fc9414ff4da72f..4d529eab1a94baf82d5b5edfac087330ffaf147e 100644
--- a/dumux/common/geometry/intersectingentities.hh
+++ b/dumux/common/geometry/intersectingentities.hh
@@ -18,378 +18,15 @@
* \file
* \ingroup Geometry
* \brief Algorithms that finds which geometric entites intersect
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_INTERSECTING_ENTITIES_HH
-#define DUMUX_INTERSECTING_ENTITIES_HH
+#ifndef DUMUX_COMMON_GEOMETRY_INTERSECTING_ENTITIES_HH
+#define DUMUX_COMMON_GEOMETRY_INTERSECTING_ENTITIES_HH
-#include <cmath>
-#include <type_traits>
-#include <vector>
-#include <algorithm>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/intersectingentities.hh"
-#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/triangulation.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief An intersection object resulting from the intersection of two primitives in an entity set
- * \note this is used as return type for some of the intersectingEntities overloads below
- */
-template<int dimworld, class CoordTypeA, class CoordTypeB = CoordTypeA>
-class IntersectionInfo
-{
-public:
- using ctype = typename Dune::PromotionTraits<CoordTypeA, CoordTypeB>::PromotedType;
- static constexpr int dimensionworld = dimworld;
- using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
-
- template<class Corners>
- explicit IntersectionInfo(std::size_t a, std::size_t b, Corners&& c)
- : a_(a)
- , b_(b)
- , corners_(c.begin(), c.end())
- {}
-
- //! 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
-};
-
-/*!
- * \ingroup Geometry
- * \brief 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,
- bool isCartesianGrid = false)
-{
- // Call the recursive find function to find candidates
- std::vector<std::size_t> entities;
- intersectingEntities(point, tree, tree.numBoundingBoxes() - 1, entities, isCartesianGrid);
- return entities;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute intersections with point for all nodes of the bounding box tree recursively
- */
-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,
- bool isCartesianGrid = false)
-{
- // 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;
- // for structured cube grids skip the primitive test
- if (isCartesianGrid)
- entities.push_back(entityIdx);
- else
- {
- 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, isCartesianGrid);
- intersectingEntities(point, tree, bBox.child1, entities, isCartesianGrid);
- }
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute all intersections between a geometry and a bounding box tree
- */
-template<class Geometry, class EntitySet>
-inline std::vector<IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>>
-intersectingEntities(const Geometry& geometry,
- const BoundingBoxTree<EntitySet>& tree)
-{
- // check if the world dimensions match
- static_assert(int(Geometry::coorddimension) == int(EntitySet::dimensionworld),
- "Can only intersect geometry and bounding box tree of same world dimension");
-
- // Create data structure for return type
- std::vector<IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>> intersections;
- using ctype = typename IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>::ctype;
- static constexpr int dimworld = Geometry::coorddimension;
-
- // compute the bounding box of the given geometry
- std::array<ctype, 2*Geometry::coorddimension> bBox;
- ctype* xMin = bBox.data(); ctype* xMax = xMin + Geometry::coorddimension;
-
- // 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]);
- }
- }
-
- // Call the recursive find function to find candidates
- intersectingEntities(geometry, tree,
- bBox, tree.numBoundingBoxes() - 1,
- intersections);
-
- return intersections;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute intersections with point for all nodes of the bounding box tree recursively
- */
-template<class Geometry, class EntitySet>
-void intersectingEntities(const Geometry& geometry,
- const BoundingBoxTree<EntitySet>& tree,
- const std::array<typename Geometry::ctype, 2*Geometry::coorddimension>& bBox,
- std::size_t nodeIdx,
- std::vector<IntersectionInfo<Geometry::coorddimension,
- typename Geometry::ctype,
- typename EntitySet::ctype>>& intersections)
-{
- // if the two bounding boxes don't intersect we can stop searching
- static constexpr int dimworld = Geometry::coorddimension;
- if (!intersectsBoundingBoxBoundingBox<dimworld>(bBox.data(), tree.getBoundingBoxCoordinates(nodeIdx)))
- return;
-
- // get node info for current bounding box node
- const auto& bBoxNode = tree.getBoundingBoxNode(nodeIdx);
-
- // if the box is a leaf do the primitive test.
- if (tree.isLeaf(bBoxNode, nodeIdx))
- {
- // eIdxA is always 0 since we intersect with exactly one geometry
- const auto eIdxA = 0;
- const auto eIdxB = bBoxNode.child1;
-
- const auto geometryTree = tree.entitySet().entity(eIdxB).geometry();
- using GeometryTree = std::decay_t<decltype(geometryTree)>;
- using Policy = IntersectionPolicy::DefaultPolicy<Geometry, GeometryTree>;
- using IntersectionAlgorithm = GeometryIntersection<Geometry, GeometryTree, Policy>;
- using Intersection = typename IntersectionAlgorithm::Intersection;
- Intersection intersection;
-
- if (IntersectionAlgorithm::intersection(geometry, geometryTree, intersection))
- {
- static constexpr int dimIntersection = Policy::dimIntersection;
- if (dimIntersection >= 2)
- {
- const auto triangulation = triangulate<dimIntersection, dimworld>(intersection);
- for (unsigned int i = 0; i < triangulation.size(); ++i)
- intersections.emplace_back(eIdxA, eIdxB, std::move(triangulation[i]));
- }
- else
- intersections.emplace_back(eIdxA, eIdxB, intersection);
- }
- }
-
- // No leaf. Check both children nodes.
- else
- {
- intersectingEntities(geometry, tree, bBox, bBoxNode.child0, intersections);
- intersectingEntities(geometry, tree, bBox, bBoxNode.child1, intersections);
- }
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute all intersections between two bounding box trees
- */
-template<class EntitySet0, class EntitySet1>
-inline std::vector<IntersectionInfo<EntitySet0::dimensionworld, typename EntitySet0::ctype, typename EntitySet1::ctype>>
-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<IntersectionInfo<EntitySet0::dimensionworld, typename EntitySet0::ctype, typename EntitySet1::ctype>> intersections;
-
- // Call the recursive find function to find candidates
- intersectingEntities(treeA, treeB,
- treeA.numBoundingBoxes() - 1,
- treeB.numBoundingBoxes() - 1,
- intersections);
-
- return intersections;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute all intersections between two all 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<IntersectionInfo<EntitySet0::dimensionworld,
- typename EntitySet0::ctype,
- typename EntitySet1::ctype>>& 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 GeometryA = std::decay_t<decltype(geometryA)>;
- using GeometryB = std::decay_t<decltype(geometryB)>;
- using Policy = IntersectionPolicy::DefaultPolicy<GeometryA, GeometryB>;
- using IntersectionAlgorithm = GeometryIntersection<GeometryA, GeometryB, Policy>;
- using Intersection = typename IntersectionAlgorithm::Intersection;
- Intersection intersection;
-
- if (IntersectionAlgorithm::intersection(geometryA, geometryB, intersection))
- {
- static constexpr int dimIntersection = Policy::dimIntersection;
-
- if (dimIntersection >= 2)
- {
- const auto triangulation = triangulate<dimIntersection, dimworld>(intersection);
- for (unsigned int i = 0; i < triangulation.size(); ++i)
- intersections.emplace_back(eIdxA, eIdxB, std::move(triangulation[i]));
- }
- else
- intersections.emplace_back(eIdxA, eIdxB, 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);
- }
-}
-
-/*!
- * \ingroup Geometry
- * \brief Compute the index of the intersecting element of a Cartesian grid with a point
- * The grid is given by the lower left corner (min), the upper right corner (max)
- * and the number of cells in each direction (cells).
- * \note If there are several options the lowest matching cell index will be returned
- * \note Returns the index i + I*j + I*J*k for the intersecting element i,j,k of a grid with I,J,K cells
- */
-template<class ctype, int dimworld>
-inline std::size_t intersectingEntityCartesianGrid(const Dune::FieldVector<ctype, dimworld>& point,
- const Dune::FieldVector<ctype, dimworld>& min,
- const Dune::FieldVector<ctype, dimworld>& max,
- const std::array<int, std::size_t(dimworld)>& cells)
-{
- std::size_t index = 0;
- for (int i = 0; i < dimworld; ++i)
- {
- using std::clamp; using std::floor;
- ctype dimOffset = clamp<ctype>(floor((point[i]-min[i])*cells[i]/(max[i]-min[i])), 0.0, cells[i]-1);
- for (int j = 0; j < i; ++j)
- dimOffset *= cells[j];
- index += static_cast<std::size_t>(dimOffset);
- }
- return index;
-}
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/intersectingentities.hh>
#endif
diff --git a/dumux/common/geometry/intersectionentityset.hh b/dumux/common/geometry/intersectionentityset.hh
index 7fad35dd08ab836297d62ac6a3880c064f38dfaa..6aa4b6f66f369b33af20da9948e86065b09e860e 100644
--- a/dumux/common/geometry/intersectionentityset.hh
+++ b/dumux/common/geometry/intersectionentityset.hh
@@ -20,271 +20,16 @@
* \file
* \ingroup Geometry
* \brief A class representing the intersection entites two geometric entity sets
+ * DEPRECATED will be removed once this header is removed
*/
#ifndef DUMUX_COMMON_GEOMETRY_INTERSECTION_ENTITY_SET_HH
#define DUMUX_COMMON_GEOMETRY_INTERSECTION_ENTITY_SET_HH
-#include <algorithm>
-#include <cassert>
-#include <iostream>
-#include <memory>
-#include <type_traits>
-#include <utility>
-#include <vector>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/intersectionentityset.hh"
-#include <dune/common/indices.hh>
-#include <dune/common/timer.hh>
-#include <dune/common/iteratorrange.hh>
-#include <dune/common/promotiontraits.hh>
-#include <dune/common/reservedvector.hh>
-#include <dune/geometry/affinegeometry.hh>
-#include <dune/geometry/type.hh>
-
-#include <dumux/common/geometry/boundingboxtree.hh>
-#include <dumux/common/geometry/intersectingentities.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief A class representing the intersection entites two geometric entity sets
- */
-template<class DomainEntitySet, class TargetEntitySet>
-class IntersectionEntitySet
-{
- using DomainTree = BoundingBoxTree<DomainEntitySet>;
- using TargetTree = BoundingBoxTree<TargetEntitySet>;
-
- using ctype = typename Dune::PromotionTraits<typename DomainEntitySet::ctype, typename TargetEntitySet::ctype>::PromotedType;
-
- static constexpr int dimWorld = DomainEntitySet::dimensionworld;
- static_assert(dimWorld == int(TargetEntitySet::dimensionworld), "Entity sets must have the same world dimension");
-
- using GlobalPosition = Dune::FieldVector<ctype, dimWorld>;
-
- static constexpr int dimDomain = DomainEntitySet::Entity::Geometry::mydimension;
- static constexpr int dimTarget = TargetEntitySet::Entity::Geometry::mydimension;
- static constexpr bool isMixedDimensional = dimDomain != dimTarget;
-
- /*!
- * \brief A class representing an intersection entity
- */
- class IntersectionEntity
- {
- static constexpr int dimIs = std::min(dimDomain, dimTarget);
- using Geometry = Dune::AffineGeometry<ctype, dimIs, dimWorld>; // geometries are always simplices
-
- // we can only have multiple neighbors in the mixeddimensional case and then only for the side with the largest dimension
- using IndexStorage = std::pair<std::conditional_t<dimDomain <= dimTarget, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>,
- std::conditional_t<dimTarget <= dimDomain, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>>;
-
- static constexpr auto domainIdx = Dune::index_constant<0>{};
- static constexpr auto targetIdx = Dune::index_constant<1>{};
-
- public:
- IntersectionEntity(const DomainTree& domainTree, const TargetTree& targetTree)
- : domainTree_(domainTree)
- , targetTree_(targetTree)
- {}
-
- //! set the intersection geometry corners
- void setCorners(const std::vector<GlobalPosition>& corners)
- {
- corners_ = corners;
- assert(corners.size() == dimIs + 1); // Only simplex intersections are allowed
- }
-
- //! add a pair of neighbor elements
- void addNeighbors(std::size_t domain, std::size_t target)
- {
- if (numDomainNeighbors() == 0 && numTargetNeighbors() == 0)
- {
- std::get<domainIdx>(neighbors_).push_back(domain);
- std::get<targetIdx>(neighbors_).push_back(target);
- }
- else if (dimDomain > dimTarget)
- std::get<domainIdx>(neighbors_).push_back(domain);
-
- else if (dimTarget > dimDomain)
- std::get<targetIdx>(neighbors_).push_back(target);
-
- else
- DUNE_THROW(Dune::InvalidStateException, "Cannot add more than one neighbor per side for equidimensional intersection!");
- }
-
- //! get the intersection geometry
- Geometry geometry() const
- { return Geometry(Dune::GeometryTypes::simplex(dimIs), corners_); }
-
- //! get the number of domain neighbors of this intersection
- std::size_t numDomainNeighbors() const
- { return std::get<domainIdx>(neighbors_).size(); }
-
- //! get the number of target neighbors of this intersection
- std::size_t numTargetNeighbors() const
- { return std::get<targetIdx>(neighbors_).size(); }
-
- //! get the nth domain neighbor entity
- typename DomainEntitySet::Entity domainEntity(unsigned int n = 0) const
- { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); }
-
- //! get the nth target neighbor entity
- typename TargetEntitySet::Entity targetEntity(unsigned int n = 0) const
- { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
-
- private:
- IndexStorage neighbors_;
- std::vector<GlobalPosition> corners_;
-
- const DomainTree& domainTree_;
- const TargetTree& targetTree_;
- };
-
- using Intersections = std::vector<IntersectionEntity>;
-
-public:
- //! make intersection entity type available
- using Entity = IntersectionEntity;
- //! make entity iterator type available
- using EntityIterator = typename Intersections::const_iterator;
-
- /*!
- * \brief Default constructor
- */
- IntersectionEntitySet() = default;
-
- /*!
- * \brief Build intersections
- */
- void build(std::shared_ptr<const DomainEntitySet> domainSet, std::shared_ptr<const TargetEntitySet> targetSet)
- {
- domainTree_ = std::make_shared<DomainTree>(domainSet);
- targetTree_ = std::make_shared<TargetTree>(targetSet);
- build(*domainTree_, *targetTree_);
- }
-
- /*!
- * \brief Build intersections
- */
- void build(std::shared_ptr<const DomainTree> domainTree, std::shared_ptr<const TargetTree> targetTree)
- {
- // make sure the tree don't get out of scope
- domainTree_ = domainTree;
- targetTree_ = targetTree;
- build(*domainTree_, *targetTree_);
- }
-
- /*!
- * \brief Build intersections
- * \note If you call this, make sure the bounding box tree stays alive for the life-time of this object
- */
- void build(const DomainTree& domainTree, const TargetTree& targetTree)
- {
- Dune::Timer timer;
-
- // compute raw intersections
- const auto rawIntersections = intersectingEntities(domainTree, targetTree);
-
- // create a map to check whether intersection geometries were already inserted
- // Note that this can only occur if the grids have different dimensionality.
- // If this is the case, we keep track of the intersections using the indices of the lower-
- // dimensional entities which is identical for all geometrically identical intersections.
- std::vector<std::vector<std::vector<GlobalPosition>>> intersectionMap;
- std::vector<std::vector<std::size_t>> intersectionIndex;
- if constexpr (isMixedDimensional)
- {
- const auto numLowDimEntities = dimTarget < dimDomain ? targetTree.entitySet().size()
- : domainTree.entitySet().size();
- intersectionMap.resize(numLowDimEntities);
- intersectionIndex.resize(numLowDimEntities);
- }
-
- // reserve memory for storing the intersections. In case of grids of
- // different dimensionality this might be an overestimate. We get rid
- // of the overhead memory at the end of this function though.
- intersections_.clear();
- intersections_.reserve(rawIntersections.size());
-
- for (const auto& rawIntersection : rawIntersections)
- {
- bool add = true;
-
- // Check if intersection was already inserted.
- // In this case we only add new neighbor information as the geometry is identical.
- if constexpr (isMixedDimensional)
- {
- const auto lowDimNeighborIdx = getLowDimNeighborIdx_(rawIntersection);
- for (int i = 0; i < intersectionMap[lowDimNeighborIdx].size(); ++i)
- {
- if (rawIntersection.cornersMatch(intersectionMap[lowDimNeighborIdx][i]))
- {
- add = false;
- // only add the pair of neighbors using the insertionIndex
- auto idx = intersectionIndex[lowDimNeighborIdx][i];
- intersections_[idx].addNeighbors(rawIntersection.first(), rawIntersection.second());
- break;
- }
- }
- }
-
- if(add)
- {
- // maybe add to the map
- if constexpr (isMixedDimensional)
- {
- intersectionMap[getLowDimNeighborIdx_(rawIntersection)].push_back(rawIntersection.corners());
- intersectionIndex[getLowDimNeighborIdx_(rawIntersection)].push_back(intersections_.size());
- }
-
- // add new intersection and add the neighbors
- intersections_.emplace_back(domainTree, targetTree);
- intersections_.back().setCorners(rawIntersection.corners());
- intersections_.back().addNeighbors(rawIntersection.first(), rawIntersection.second());
- }
- }
-
- intersections_.shrink_to_fit();
- std::cout << "Computed " << size() << " intersection entities in " << timer.elapsed() << std::endl;
- }
-
- //! return begin iterator to intersection container
- typename Intersections::const_iterator ibegin() const
- { return intersections_.begin(); }
-
- //! return end iterator to intersection container
- typename Intersections::const_iterator iend() const
- { return intersections_.end(); }
-
- //! the number of intersections
- std::size_t size() const
- { return intersections_.size(); }
-
- /*!
- * \brief Range generator to iterate with range-based for loops over all intersections
- * as follows: for (const auto& is : intersections(glue)) { ... }
- * \note free function
- */
- friend Dune::IteratorRange<EntityIterator> intersections(const IntersectionEntitySet& set)
- { return {set.ibegin(), set.iend()}; }
-
-private:
- template<class RawIntersection,
- bool enable = isMixedDimensional, std::enable_if_t<enable, int> = 0>
- auto getLowDimNeighborIdx_(const RawIntersection& is)
- {
- if constexpr (dimTarget < dimDomain)
- return is.second();
- else
- return is.first();
- }
-
- Intersections intersections_;
-
- std::shared_ptr<const DomainTree> domainTree_;
- std::shared_ptr<const TargetTree> targetTree_;
-};
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/intersectionentityset.hh>
#endif
diff --git a/dumux/common/geometry/intersectspointgeometry.hh b/dumux/common/geometry/intersectspointgeometry.hh
index 68f0c8b115da3f04f6722e48973bca68d2efee91..909c6d1b58ae08ed67f3f014b9be2cd50b10b48c 100644
--- a/dumux/common/geometry/intersectspointgeometry.hh
+++ b/dumux/common/geometry/intersectspointgeometry.hh
@@ -18,104 +18,15 @@
* \file
* \ingroup Geometry
* \brief Detect if a point intersects a geometry
+ * DEPRECATED will be removed once this header is removed
*/
#ifndef DUMUX_INTERSECTS_POINT_GEOMETRY_HH
#define DUMUX_INTERSECTS_POINT_GEOMETRY_HH
-#include <cmath>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/intersectspointgeometry.hh"
-#include "intersectspointsimplex.hh"
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief 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.");
-}
-
-/*!
- * \ingroup Geometry
- * \brief 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.");
-}
-
-/*!
- * \ingroup Geometry
- * \brief 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
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/intersectspointgeometry.hh>
#endif
diff --git a/dumux/common/geometry/intersectspointsimplex.hh b/dumux/common/geometry/intersectspointsimplex.hh
index d3eb254c37729eee6a5b52b37586fb987f05ff85..cb8639f536f0dc8617ce786f92dfec7ae321a0f8 100644
--- a/dumux/common/geometry/intersectspointsimplex.hh
+++ b/dumux/common/geometry/intersectspointsimplex.hh
@@ -18,208 +18,14 @@
* \file
* \ingroup Geometry
* \brief Detect if a point intersects a simplex (including boundary)
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_INTERSECTS_POINT_SIMPLEX_HH
-#define DUMUX_INTERSECTS_POINT_SIMPLEX_HH
+#ifndef DUMUX_COMMON_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
+#define DUMUX_COMMON_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
-#include <cmath>
-#include <dune/common/fvector.hh>
-#include <dumux/common/math.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief 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;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Find out whether a point is inside a triangle (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)
-{
- // adapted from the algorithm from from "Real-Time Collision Detection" by Christer Ericson,
- // published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.4.2)
- constexpr ctype eps_ = 1.0e-7;
-
- // compute the normal of the triangle
- const auto v1 = p0 - p2;
- auto n = crossProduct(v1, p1 - p0);
- const ctype nnorm = n.two_norm();
- const ctype eps4 = eps_*nnorm*nnorm; // compute an epsilon for later
- n /= nnorm; // normalize
-
- // first check if we are in the plane of the triangle
- // if not we can return early
- using std::abs;
- auto x = p0 - point;
- x /= x.two_norm(); // normalize
-
- if (abs(x*n) > eps_)
- return false;
-
- // translate the triangle so that 'point' is the origin
- const auto a = p0 - point;
- const auto b = p1 - point;
- const auto c = p2 - point;
-
- // compute the normal vectors for triangles P->A->B and P->B->C
- const auto u = crossProduct(b, c);
- const auto v = crossProduct(c, a);
-
- // they have to point in the same direction or be orthogonal
- if (u*v < 0.0 - eps4)
- return false;
-
- // compute the normal vector for triangle P->C->A
- const auto w = crossProduct(a, b);
-
- // it also has to point in the same direction or be orthogonal
- if (u*w < 0.0 - eps4)
- return false;
-
- // check if point is on the line of one of the edges
- if (u.two_norm2() < eps4)
- return b*c < 0.0 + eps_*nnorm;
- if (v.two_norm2() < eps4)
- return a*c < 0.0 + eps_*nnorm;
-
- // now the point must be in the triangle (or on the faces)
- return true;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Find out whether a point is inside a triangle (p0, p1, p2, p3) (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);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Find out whether a point is inside a interval (p0, p1, p2, p3) (dimworld is 3)
- * \note We assume the given interval has non-zero length and use it to scale the epsilon
- */
-template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3 || 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;
-
- const ctype v1norm = v1.two_norm();
- const ctype v2norm = v2.two_norm();
-
- // check if point and p0 are the same
- if (v2norm < v1norm*eps_)
- return true;
-
- return (v1.dot(v2) > v1norm*v2norm*(1.0 - eps_) && v2norm < v1norm*(1.0 + eps_));
-}
-
-/*!
- * \ingroup Geometry
- * \brief Find out whether a point is inside a interval (p0, p1, p2, p3) (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
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/intersectspointsimplex.hh"
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/intersectspointsimplex.hh>
#endif
diff --git a/dumux/common/geometry/makegeometry.hh b/dumux/common/geometry/makegeometry.hh
index 68a057e53204c65b3e82614dd0c741245f7b7ca1..2590962b50b9e04ff507b2fe51978f0597e12a08 100644
--- a/dumux/common/geometry/makegeometry.hh
+++ b/dumux/common/geometry/makegeometry.hh
@@ -18,193 +18,15 @@
* \file
* \ingroup Geometry
* \brief Create Dune geometries from user-specified points
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_MAKE_GEOMETRY_HH
-#define DUMUX_MAKE_GEOMETRY_HH
+#ifndef DUMUX_COMMON_GEOMETRY_MAKE_GEOMETRY_HH
+#define DUMUX_COMMON_GEOMETRY_MAKE_GEOMETRY_HH
-#include <vector>
-#include <array>
-#include <limits>
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/geometry/multilineargeometry.hh>
-#include <dumux/common/math.hh>
-#include <dumux/common/geometry/grahamconvexhull.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/makegeometry.hh"
-namespace Dumux {
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/makegeometry.hh>
-/*!
- * \ingroup Geometry
- * \brief Checks if four points lie within the same plane
- */
-template<class CoordScalar>
-bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points, const CoordScalar scale)
-{
- if (points.size() != 4)
- DUNE_THROW(Dune::InvalidStateException, "Check only works for 4 points!");
-
- // (see "Real-Time Collision Detection" by Christer Ericson)
- Dune::FieldMatrix<CoordScalar, 4, 4> M;
- for (int i = 0; i < 3; ++i )
- M[i] = {points[0][i], points[1][i], points[2][i], points[3][i]};
- M[3] = {1.0*scale, 1.0*scale, 1.0*scale, 1.0*scale};
-
- using std::abs;
- return abs(M.determinant()) < 1.5e-7*scale*scale*scale*scale;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Checks if four points lie within the same plane.
- */
-template<class CoordScalar>
-bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
-{
- Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
- Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
- for (const auto& p : points)
- {
- for (int i=0; i<3; i++)
- {
- using std::min;
- using std::max;
- bBoxMin[i] = min(bBoxMin[i], p[i]);
- bBoxMax[i] = max(bBoxMax[i], p[i]);
- }
- }
-
- const auto size = (bBoxMax - bBoxMin).two_norm();
-
- return pointsAreCoplanar(points, size);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Returns a vector of points following the dune ordering.
- * Convenience method that creates a temporary object in case no array of orientations is desired.
- *
- * \param points The user-specified vector of points (potentially in wrong order).
- */
-template<class CoordScalar>
-std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
-{
- std::array<int, 4> tmp;
- return getReorderedPoints(points, tmp);
-}
-
-/*!
- * \ingroup Geometry
- * \brief Returns a vector of points following the dune ordering.
- *
- * \param points The user-specified vector of points (potentially in wrong order).
- * \param orientations An array of orientations that can be useful for further processing.
- */
-template<class CoordScalar>
-std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points,
- std::array<int, 4>& orientations)
-{
- if(points.size() == 4)
- {
- auto& p0 = points[0];
- auto& p1 = points[1];
- auto& p2 = points[2];
- auto& p3 = points[3];
-
- // check if the points define a proper quadrilateral
- const auto normal = crossProduct((p1 - p0), (p2 - p0));
-
- orientations = { getOrientation(p0, p3, p2, normal),
- getOrientation(p0, p3, p1, normal),
- getOrientation(p2, p1, p0, normal),
- getOrientation(p2, p1, p3, normal) };
-
-
- // check if the points follow the dune ordering (see http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf)
- const bool diagonalsIntersect = (orientations[0] != orientations[1]) && (orientations[2] != orientations[3]);
-
- // the points conform with the dune ordering
- if(diagonalsIntersect)
- return points;
-
- // the points do not conform with the dune ordering, re-order
- using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
- if(!diagonalsIntersect && orientations[0] == 1)
- return std::vector<GlobalPosition>{p1, p0, p2, p3};
- else if(!diagonalsIntersect && orientations[0] == -1)
- return std::vector<GlobalPosition>{p3, p1, p0, p2};
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not reorder points");
- }
- else
- DUNE_THROW(Dune::NotImplemented, "Reorder for " << points.size() << " points.");
-}
-
-/*!
- * \ingroup Geometry
- * \brief Creates a dune quadrilateral geometry given 4 corner points.
- *
- * \tparam CoordScalar The CoordScalar type.
- * \tparam enableSanityCheck Turn on/off sanity check and reordering of points
- * \param points The user-specified vector of points (potentially in wrong order).
- */
-template<class CoordScalar, bool enableSanityCheck = true>
-auto makeDuneQuadrilaterial(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
-{
- if (points.size() != 4)
- DUNE_THROW(Dune::InvalidStateException, "A quadrilateral needs 4 corner points!");
-
- using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
- static constexpr auto coordDim = GlobalPosition::dimension;
- static constexpr auto dim = coordDim-1;
- using GeometryType = Dune::MultiLinearGeometry<CoordScalar, dim, coordDim>;
-
- // if no sanity check if desired, use the given points directly to construct the geometry
- if (!enableSanityCheck)
- return GeometryType(Dune::GeometryTypes::quadrilateral, points);
-
- // compute size
- Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
- Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
- for (const auto& p : points)
- {
- for (int i = 0; i < 3; i++)
- {
- using std::min;
- using std::max;
- bBoxMin[i] = min(bBoxMin[i], p[i]);
- bBoxMax[i] = max(bBoxMax[i], p[i]);
- }
- }
-
- const auto size = (bBoxMax - bBoxMin).two_norm();
-
- // otherwise, perform a number of checks and corrections
- if (!pointsAreCoplanar(points, size))
- DUNE_THROW(Dune::InvalidStateException, "Points do not lie within a plane");
-
- auto corners = grahamConvexHull<2>(points);
- if (corners.size() != 4)
- DUNE_THROW(Dune::InvalidStateException, "Points do not span a strictly convex polygon!");
-
- // make sure points conform with dune ordering
- std::array<int, 4> orientations;
- corners = getReorderedPoints(corners, orientations);
-
- if (std::any_of(orientations.begin(), orientations.end(), [](auto i){ return i == 0; }))
- DUNE_THROW(Dune::InvalidStateException, "More than two points lie on the same line.");
-
- const auto quadrilateral = GeometryType(Dune::GeometryTypes::quadrilateral, corners);
-
- const auto eps = 1e-7;
- if (quadrilateral.volume() < eps*size*size)
- DUNE_THROW(Dune::InvalidStateException, "Something went wrong, geometry has area of zero");
-
- return quadrilateral;
-}
-
-
-
-} // end namespace Dumux
-
-#endif // DUMUX_MAKE_GEOMETRY_HH
+#endif
diff --git a/dumux/common/geometry/normal.hh b/dumux/common/geometry/normal.hh
index a6e95ecb170b8c75a6df877f8667742f85585ee7..712580c5efae33313c3765ecf60e8b3d2535186a 100644
--- a/dumux/common/geometry/normal.hh
+++ b/dumux/common/geometry/normal.hh
@@ -18,60 +18,15 @@
* \file
* \ingroup Geometry
* \brief Helper functions for normals
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_GEOMETRY_NORMAL_HH
-#define DUMUX_GEOMETRY_NORMAL_HH
+#ifndef DUMUX_COMMON_GEOMETRY_NORMAL_HH
+#define DUMUX_COMMON_GEOMETRY_NORMAL_HH
-#include <algorithm>
-#include <dune/common/float_cmp.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/normal.hh"
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief Create a vector normal to the given one (v is expected to be non-zero)
- * \note This returns some orthogonal vector with arbitrary length
- */
-template<class Vector>
-inline Vector normal(const Vector& v)
-{
- static_assert(Vector::size() > 1, "normal expects a coordinate dimension > 1");
-
- if constexpr (Vector::size() == 2)
- return Vector({-v[1], v[0]});
-
- const auto it = std::find_if(v.begin(), v.end(), [](const auto& x) { return Dune::FloatCmp::ne(x, 0.0); });
- const auto index = std::distance(v.begin(), it);
- if (index != Vector::size()-1)
- {
- Vector normal(0.0);
- normal[index] = -v[index+1];
- normal[index+1] = v[index];
- return normal;
- }
- else
- {
- Vector normal(0.0);
- normal[index-1] = -v[index];
- normal[index] = v[index-1];
- return normal;
-
- }
-}
-
-/*!
- * \ingroup Geometry
- * \brief Create a vector normal to the given one (v is expected to be non-zero)
- * \note This returns some orthogonal vector with unit length
- */
-template<class Vector>
-inline Vector unitNormal(const Vector& v)
-{
- auto normal = Dumux::normal(v);
- normal /= normal.two_norm();
- return normal;
-}
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/normal.hh>
#endif
diff --git a/dumux/common/geometry/refinementquadraturerule.hh b/dumux/common/geometry/refinementquadraturerule.hh
index d2996c623ef0c64b8b7c2ab37c56c12ca80fcc71..0778c474a0d5a9ee06bf920204f048d19ffa15b5 100644
--- a/dumux/common/geometry/refinementquadraturerule.hh
+++ b/dumux/common/geometry/refinementquadraturerule.hh
@@ -21,74 +21,16 @@
* \ingroup Geometry
* \author Timo Koch
* \brief A quadrature based on refinement
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_COMMON_REFINEMENT_QUADRATURERULE_HH
-#define DUMUX_COMMON_REFINEMENT_QUADRATURERULE_HH
+#ifndef DUMUX_COMMON_GEOMETRY_REFINEMENT_QUADRATURERULE_HH
+#define DUMUX_COMMON_GEOMETRY_REFINEMENT_QUADRATURERULE_HH
-#include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/virtualrefinement.hh>
-#include <dumux/common/math.hh>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/refinementquadraturerule.hh"
-namespace Dumux {
-
-/*!
- * \ingroup Geometry
- * \brief A "quadrature" based on virtual refinement
- */
-template<typename ct, int mydim>
-class RefinementQuadratureRule final : public Dune::QuadratureRule<ct, mydim>
-{
-public:
- //! The space dimension
- static constexpr int dim = mydim;
-
- //! The highest quadrature order available
- static constexpr int highest_order = 1;
-
- ~RefinementQuadratureRule() final {}
-
- RefinementQuadratureRule(Dune::GeometryType type, int levels)
- : Dune::QuadratureRule<ct, dim>(type, highest_order)
- {
- auto& refinement = Dune::buildRefinement<dim, ct>(type, type);
- const auto tag = Dune::refinementLevels(levels);
-
- // get the vertices
- const auto numVertices = refinement.nVertices(tag);
- std::vector<Dune::FieldVector<ct, dim>> points; points.reserve(numVertices);
- auto vIt = refinement.vBegin(tag);
- const auto vItEnd = refinement.vEnd(tag);
- for (; vIt != vItEnd; ++vIt)
- points.emplace_back(vIt.coords());
-
- // go over all elements and get center and volume
- const auto numElements = refinement.nElements(tag);
- this->reserve(numElements);
- auto eIt = refinement.eBegin(tag);
- const auto eItEnd = refinement.eEnd(tag);
- for (; eIt != eItEnd; ++eIt)
- {
- // quadrature point in the centroid of the sub-element and weight is its volume
- const auto weight = computeVolume_(type, points, eIt.vertexIndices());
- this->emplace_back(eIt.coords(), weight);
- }
- }
-
-private:
- template<class VertexIndices>
- ct computeVolume_(Dune::GeometryType t, const std::vector<Dune::FieldVector<ct, dim>>& points, const VertexIndices& indices) const
- {
- if (t != Dune::GeometryTypes::simplex(2))
- DUNE_THROW(Dune::NotImplemented, "Only implemented for 2d simplices");
-
- const auto ab = points[indices[1]] - points[indices[0]];
- const auto ac = points[indices[2]] - points[indices[0]];
- using std::abs;
- return abs(crossProduct(ab, ac))*0.5;
- }
-};
-
-} // end namespace Dumux
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/refinementquadraturerule.hh>
#endif
diff --git a/dumux/common/geometry/triangulation.hh b/dumux/common/geometry/triangulation.hh
index 2b04870cd6e7357f7d9ebc18a1fd2323caac982a..1cba28a5832a197fe6d09cbce8d3ec216fe2f7be 100644
--- a/dumux/common/geometry/triangulation.hh
+++ b/dumux/common/geometry/triangulation.hh
@@ -18,117 +18,15 @@
* \file
* \ingroup Geometry
* \brief Functionality to triangulate point clouds
+ * DEPRECATED will be removed once this header is removed
*/
-#ifndef DUMUX_TRIANGULATION_HH
-#define DUMUX_TRIANGULATION_HH
+#ifndef DUMUX_COMMON_GEOMETRY_TRIANGULATION_HH
+#define DUMUX_COMMON_GEOMETRY_TRIANGULATION_HH
-#include <vector>
-#include <array>
-#include <algorithm>
-#include <type_traits>
+#warning "This header is deprecated and will be removed after release 3.3. Please use dumux/geometry/triangulation.hh"
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
+// This header, and all other geometry headers have been moved to their own folder.
+// Please use the geometry headers in dumux/geometry/, as this will be removed after release 3.3.
+#include <dumux/geometry/triangulation.hh>
-#include <dumux/common/math.hh>
-
-namespace Dumux {
-namespace TriangulationPolicy {
-
-//! Policy that expects a point cloud that represents a convex
-//! hull. Inserts the mid point and connects it to the points of the hull.
-struct MidPointPolicy {};
-
-//! Delaunay-type triangulations.
-struct DelaunayPolicy {};
-
-template<int dim, int dimWorld>
-using DefaultPolicy = std::conditional_t< dim >= 2, MidPointPolicy, DelaunayPolicy >;
-
-} // end namespace TriangulationPolicy
-
-//! The data type to store triangulations
-template<int dim, int dimWorld, class ctype>
-using Triangulation = std::vector< std::array<Dune::FieldVector<ctype, dimWorld>, dim+1> >;
-
-/*!
- * \ingroup Geometry
- * \brief Triangulate area given points of a convex hull
- * \tparam dim Specifies the dimension of the resulting triangulation (2 or 3)
- * \tparam dimWorld The dimension of the coordinate space
- * \tparam Policy Specifies the algorithm to be used for triangulation
- *
- * \note this specialization is for 2d triangulations using mid point policy
- * \note Assumes all points of the convex hull are coplanar
- * \note This inserts a mid point and connects all corners with that point to triangles
- */
-template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
- class RandomAccessContainer,
- std::enable_if_t< std::is_same<Policy, TriangulationPolicy::MidPointPolicy>::value
- && dim == 2, int> = 0 >
-inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type>
-triangulate(const RandomAccessContainer& convexHull)
-{
- using ctype = typename RandomAccessContainer::value_type::value_type;
- using Point = Dune::FieldVector<ctype, dimWorld>;
- using Triangle = std::array<Point, 3>;
-
- static_assert(std::is_same<typename RandomAccessContainer::value_type, Point>::value,
- "Triangulation expects Dune::FieldVector as point type");
-
- if (convexHull.size() < 3)
- DUNE_THROW(Dune::InvalidStateException, "Try to triangulate point cloud with less than 3 points!");
-
- if (convexHull.size() == 3)
- return std::vector<Triangle>(1, {convexHull[0], convexHull[1], convexHull[2]});
-
- Point midPoint(0.0);
- for (const auto p : convexHull)
- midPoint += p;
- midPoint /= convexHull.size();
-
- std::vector<Triangle> triangulation;
- triangulation.reserve(convexHull.size());
-
- for (std::size_t i = 0; i < convexHull.size()-1; ++i)
- triangulation.emplace_back(Triangle{midPoint, convexHull[i], convexHull[i+1]});
-
- triangulation.emplace_back(Triangle{midPoint, convexHull[convexHull.size()-1], convexHull[0]});
-
- return triangulation;
-}
-
-/*!
- * \ingroup Geometry
- * \brief Triangulate area given points of a convex hull
- * \tparam dim Specifies the dimension of the resulting triangulation (2 or 3)
- * \tparam dimWorld The dimension of the coordinate space
- * \tparam Policy Specifies the algorithm to be used for triangulation
- *
- * \note this specialization is for 1d discretization using segments
- * \note sorts the points and connects them via segments
- */
-template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
- class RandomAccessContainer,
- std::enable_if_t< std::is_same<Policy, TriangulationPolicy::DelaunayPolicy>::value
- && dim == 1, int> = 0 >
-inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type>
-triangulate(const RandomAccessContainer& points)
-{
- using ctype = typename RandomAccessContainer::value_type::value_type;
- using Point = Dune::FieldVector<ctype, dimWorld>;
-
- static_assert(std::is_same<typename RandomAccessContainer::value_type, Point>::value,
- "Triangulation expects Dune::FieldVector as point type");
-
- if (points.size() == 2)
- return Triangulation<dim, dimWorld, ctype>({ {points[0], points[1]} });
-
- //! \todo sort points and create polyline
- assert(points.size() > 1);
- DUNE_THROW(Dune::NotImplemented, "1d triangulation for point cloud size > 2");
-}
-
-} // end namespace Dumux
-
-# endif
+#endif