diff --git a/dumux/geometry/intersectionentityset.hh b/dumux/geometry/intersectionentityset.hh
index 7152297dec2cde2d309547f02517e2f91d85bb51..6f328df4b57bef3556ccdb6b278f836f0b25720c 100644
--- a/dumux/geometry/intersectionentityset.hh
+++ b/dumux/geometry/intersectionentityset.hh
@@ -32,105 +32,123 @@
#include <dumux/geometry/boundingboxtree.hh>
#include <dumux/geometry/intersectingentities.hh>
-namespace Dumux {
+namespace Dumux::Detail::Intersections {
/*!
- * \ingroup Geometry
- * \brief A class representing the intersection entities two geometric entity sets
+ * \brief A class representing an intersection entity
*/
template<class DomainEntitySet, class TargetEntitySet>
-class IntersectionEntitySet
+class IntersectionEntity
{
- using DomainTree = BoundingBoxTree<DomainEntitySet>;
- using TargetTree = BoundingBoxTree<TargetEntitySet>;
-
+ static constexpr int dimDomain = DomainEntitySet::Entity::Geometry::mydimension;
+ static constexpr int dimTarget = TargetEntitySet::Entity::Geometry::mydimension;
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");
+ static constexpr int dimIs = std::min(dimDomain, dimTarget);
+
using GlobalPosition = Dune::FieldVector<ctype, dimWorld>;
+ using Geometry = Dune::AffineGeometry<ctype, dimIs, dimWorld>; // geometries are always simplices
- static constexpr int dimDomain = DomainEntitySet::Entity::Geometry::mydimension;
- static constexpr int dimTarget = TargetEntitySet::Entity::Geometry::mydimension;
- static constexpr bool isMixedDimensional = dimDomain != dimTarget;
+ // 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>>>;
- /*!
- * \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
+ static constexpr auto domainIdx = Dune::index_constant<0>{};
+ static constexpr auto targetIdx = Dune::index_constant<1>{};
- // 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 DomainEntitySet& domainEntitySet, const TargetEntitySet& targetEntitySet)
+ : domainEntitySet_(domainEntitySet)
+ , targetEntitySet_(targetEntitySet)
+ {}
- 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
+ }
- //! set the intersection geometry corners
- void setCorners(const std::vector<GlobalPosition>& corners)
+ //! add a pair of neighbor elements
+ void addNeighbors(std::size_t domain, std::size_t target)
+ {
+ if (numDomainNeighbors() == 0 && numTargetNeighbors() == 0)
{
- corners_ = corners;
- assert(corners.size() == dimIs + 1); // Only simplex intersections are allowed
+ 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);
- //! 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 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!");
+ }
- 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 domainEntitySet_.entity(std::get<domainIdx>(neighbors_)[n]); }
+
+ //! get the nth target neighbor entity
+ typename TargetEntitySet::Entity targetEntity(unsigned int n = 0) const
+ { return targetEntitySet_.entity(std::get<targetIdx>(neighbors_)[n]); }
- //! get the intersection geometry
- Geometry geometry() const
- { return Geometry(Dune::GeometryTypes::simplex(dimIs), corners_); }
+private:
+ IndexStorage neighbors_;
+ std::vector<GlobalPosition> corners_;
- //! get the number of domain neighbors of this intersection
- std::size_t numDomainNeighbors() const
- { return std::get<domainIdx>(neighbors_).size(); }
+ const DomainEntitySet& domainEntitySet_;
+ const TargetEntitySet& targetEntitySet_;
+};
- //! get the number of target neighbors of this intersection
- std::size_t numTargetNeighbors() const
- { return std::get<targetIdx>(neighbors_).size(); }
+} // end namespace Dumux::Detail::Intersections
- //! get the nth domain neighbor entity
- typename DomainEntitySet::Entity domainEntity(unsigned int n = 0) const
- { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); }
+namespace Dumux {
- //! get the nth target neighbor entity
- typename TargetEntitySet::Entity targetEntity(unsigned int n = 0) const
- { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
+/*!
+ * \ingroup Geometry
+ * \brief A class representing the intersection entities two geometric entity sets
+ */
+template<class DomainEntitySet, class TargetEntitySet>
+class IntersectionEntitySet
+{
+ using DomainTree = BoundingBoxTree<DomainEntitySet>;
+ using TargetTree = BoundingBoxTree<TargetEntitySet>;
- private:
- IndexStorage neighbors_;
- std::vector<GlobalPosition> corners_;
+ using ctype = typename Dune::PromotionTraits<typename DomainEntitySet::ctype, typename TargetEntitySet::ctype>::PromotedType;
- const DomainTree& domainTree_;
- const TargetTree& targetTree_;
- };
+ 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;
+
+ using IntersectionEntity = Dumux::Detail::Intersections::IntersectionEntity<DomainEntitySet, TargetEntitySet>;
using Intersections = std::vector<IntersectionEntity>;
+ using DomainGeometry = typename DomainEntitySet::Entity::Geometry;
+ using TargetGeometry = typename TargetEntitySet::Entity::Geometry;
+
public:
//! make intersection entity type available
using Entity = IntersectionEntity;
@@ -147,9 +165,26 @@ public:
*/
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_);
+ // specialized algorithm if the target geometry is only a single entity
+ if (targetSet->size() == 1)
+ {
+ domainTree_ = std::make_shared<DomainTree>(domainSet);
+ targetEntitySet_ = targetSet;
+ build_(*domainTree_, targetEntitySet_->entity(0).geometry());
+ }
+ // specialized algorithm if the domain geometry is only a single entity
+ else if (domainSet->size() == 1)
+ {
+ domainEntitySet_ = domainSet;
+ targetTree_ = std::make_shared<TargetTree>(targetSet);
+ build_(domainEntitySet_->entity(0).geometry(), *targetTree_);
+ }
+ else
+ {
+ domainTree_ = std::make_shared<DomainTree>(domainSet);
+ targetTree_ = std::make_shared<TargetTree>(targetSet);
+ build(*domainTree_, *targetTree_);
+ }
}
/*!
@@ -171,9 +206,62 @@ public:
{
Dune::Timer timer;
- // compute raw intersections
const auto rawIntersections = intersectingEntities(domainTree, targetTree);
+ build_(domainTree.entitySet(), targetTree.entitySet(), rawIntersections);
+
+ 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:
+ void build_(const DomainTree& domainTree, const TargetGeometry& targetGeometry)
+ {
+ Dune::Timer timer;
+
+ const auto rawIntersections = intersectingEntities(targetGeometry, domainTree);
+ // because in intersectingEntities geometry always comes first (domain set)
+ // but here we want the geometry to be the second (target set), we need to flip the index sets
+ // of the resulting intersections
+ const bool flipNeighborIndices = true;
+ build_(domainTree.entitySet(), *targetEntitySet_, rawIntersections, flipNeighborIndices);
+
+ std::cout << "Computed " << size() << " intersection entities in " << timer.elapsed() << std::endl;
+ }
+
+ void build_(const DomainGeometry& domainGeometry, const TargetTree& targetTree)
+ {
+ Dune::Timer timer;
+
+ const auto rawIntersections = intersectingEntities(domainGeometry, targetTree);
+ build_(*domainEntitySet_, targetTree.entitySet(), rawIntersections);
+
+ std::cout << "Computed " << size() << " intersection entities in " << timer.elapsed() << std::endl;
+ }
+ template<class RawIntersections>
+ void build_(const DomainEntitySet& domainEntitySet, const TargetEntitySet& targetEntitySet,
+ const RawIntersections& rawIntersections,
+ bool flipNeighborIndices = false)
+ {
// 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-
@@ -182,8 +270,7 @@ public:
std::vector<std::vector<std::size_t>> intersectionIndex;
if constexpr (isMixedDimensional)
{
- const auto numLowDimEntities = dimTarget < dimDomain ? targetTree.entitySet().size()
- : domainTree.entitySet().size();
+ const auto numLowDimEntities = dimTarget < dimDomain ? targetEntitySet.size() : domainEntitySet.size();
intersectionMap.resize(numLowDimEntities);
intersectionIndex.resize(numLowDimEntities);
}
@@ -202,7 +289,7 @@ public:
// In this case we only add new neighbor information as the geometry is identical.
if constexpr (isMixedDimensional)
{
- const auto lowDimNeighborIdx = getLowDimNeighborIdx_(rawIntersection);
+ const auto lowDimNeighborIdx = getLowDimNeighborIndex_(rawIntersection, flipNeighborIndices);
for (int i = 0; i < intersectionMap[lowDimNeighborIdx].size(); ++i)
{
if (rawIntersection.cornersMatch(intersectionMap[lowDimNeighborIdx][i]))
@@ -210,7 +297,7 @@ public:
add = false;
// only add the pair of neighbors using the insertionIndex
auto idx = intersectionIndex[lowDimNeighborIdx][i];
- intersections_[idx].addNeighbors(rawIntersection.first(), rawIntersection.second());
+ intersections_[idx].addNeighbors(domainNeighborIndex_(rawIntersection, flipNeighborIndices), targetNeighborIndex_(rawIntersection, flipNeighborIndices));
break;
}
}
@@ -221,56 +308,45 @@ public:
// maybe add to the map
if constexpr (isMixedDimensional)
{
- intersectionMap[getLowDimNeighborIdx_(rawIntersection)].push_back(rawIntersection.corners());
- intersectionIndex[getLowDimNeighborIdx_(rawIntersection)].push_back(intersections_.size());
+ const auto lowDimNeighborIdx = getLowDimNeighborIndex_(rawIntersection, flipNeighborIndices);
+ intersectionMap[lowDimNeighborIdx].push_back(rawIntersection.corners());
+ intersectionIndex[lowDimNeighborIdx].push_back(intersections_.size());
}
// add new intersection and add the neighbors
- intersections_.emplace_back(domainTree, targetTree);
+ intersections_.emplace_back(domainEntitySet, targetEntitySet);
intersections_.back().setCorners(rawIntersection.corners());
- intersections_.back().addNeighbors(rawIntersection.first(), rawIntersection.second());
+ intersections_.back().addNeighbors(domainNeighborIndex_(rawIntersection, flipNeighborIndices), targetNeighborIndex_(rawIntersection, flipNeighborIndices));
}
}
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)
+ auto getLowDimNeighborIndex_(const RawIntersection& is, bool flipNeighborIndices)
{
if constexpr (dimTarget < dimDomain)
- return is.second();
+ return targetNeighborIndex_(is, flipNeighborIndices);
else
- return is.first();
+ return domainNeighborIndex_(is, flipNeighborIndices);
}
+ template<class RawIntersection>
+ auto domainNeighborIndex_(const RawIntersection& is, bool flipNeighborIndices)
+ { return flipNeighborIndices ? is.second() : is.first(); }
+
+ template<class RawIntersection>
+ auto targetNeighborIndex_(const RawIntersection& is, bool flipNeighborIndices)
+ { return flipNeighborIndices ? is.first() : is.second(); }
+
Intersections intersections_;
std::shared_ptr<const DomainTree> domainTree_;
std::shared_ptr<const TargetTree> targetTree_;
+ std::shared_ptr<const DomainEntitySet> domainEntitySet_;
+ std::shared_ptr<const TargetEntitySet> targetEntitySet_;
};
} // end namespace Dumux