diff --git a/dumux/discretization/projection/projector.hh b/dumux/discretization/projection/projector.hh index 1ed7549d1575e2f7df77ebbe4ec8278c2cb82548..fe9160393b0bc205b763590ad6cd5c8c4798f9f9 100644 --- a/dumux/discretization/projection/projector.hh +++ b/dumux/discretization/projection/projector.hh @@ -334,14 +334,13 @@ auto createProjectionMatrices(const FEBasisDomain& feBasisDomain, unsigned int maxBasisOrder = 0; for (const auto& is : intersections(glue)) { - // "outside" element is of target domain // since target dim <= domain dim there is maximum one! - targetLocalView.bind( is.outside(0) ); + targetLocalView.bind( is.targetEntity(0) ); const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis(); - for (unsigned int nIdx = 0; nIdx < is.neighbor(/*targetSide*/1); ++nIdx) + for (unsigned int nIdx = 0; nIdx < is.numDomainNeighbors(); ++nIdx) { - domainLocalView.bind( is.inside(nIdx) ); + domainLocalView.bind( is.domainEntity(nIdx) ); const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis(); // keep track of maximum basis order (used in integration) @@ -370,7 +369,7 @@ auto createProjectionMatrices(const FEBasisDomain& feBasisDomain, for (const auto& is : intersections(glue)) { - const auto& targetElement = is.outside(0); + const auto& targetElement = is.targetEntity(0); const auto& targetElementGeometry = targetElement.geometry(); targetLocalView.bind( targetElement ); @@ -411,10 +410,10 @@ auto createProjectionMatrices(const FEBasisDomain& feBasisDomain, // targetElement is aligned with a facet of domainElement. In // this case make sure the basis functions are not added // multiple times! (division by numNeighbors) - const auto numNeighbors = is.neighbor(/*targetSide*/1); + const auto numNeighbors = is.numDomainNeighbors(); for (unsigned int nIdx = 0; nIdx < numNeighbors; ++nIdx) { - const auto& domainElement = is.inside(nIdx); + const auto& domainElement = is.domainEntity(nIdx); domainLocalView.bind( domainElement ); const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis(); diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh index 37a30c93e4109b25ac8f66e076e11588776c8b25..5fc41635d9b308e7bdef8993d7843226704bfc5a 100644 --- a/dumux/multidomain/embedded/couplingmanagerbase.hh +++ b/dumux/multidomain/embedded/couplingmanagerbase.hh @@ -39,7 +39,7 @@ #include <dumux/common/geometry/intersectingentities.hh> #include <dumux/discretization/method.hh> #include <dumux/multidomain/couplingmanager.hh> -#include <dumux/multidomain/embedded/mixeddimensionglue.hh> +#include <dumux/multidomain/glue.hh> #include <dumux/multidomain/embedded/pointsourcedata.hh> #include <dumux/multidomain/embedded/integrationpointsource.hh> @@ -104,7 +104,7 @@ class EmbeddedCouplingManagerBase using CouplingStencil = std::vector<std::size_t>; using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate; - using GlueType = MultiDomainGlue<GridView<lowDimIdx>, GridView<bulkIdx>, ElementMapper<lowDimIdx>, ElementMapper<bulkIdx>>; + using GlueType = MultiDomainGlue<GridView<bulkIdx>, GridView<lowDimIdx>, ElementMapper<bulkIdx>, ElementMapper<lowDimIdx>>; public: //! export traits @@ -482,7 +482,7 @@ protected: const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).fvGridGeometry(); // intersect the bounding box trees - glue_->build(lowDimFvGridGeometry.boundingBoxTree(), bulkFvGridGeometry.boundingBoxTree()); + glue_->build(bulkFvGridGeometry.boundingBoxTree(), lowDimFvGridGeometry.boundingBoxTree()); } template<class Geometry, class GlobalPosition> diff --git a/dumux/multidomain/embedded/mixeddimensionglue.hh b/dumux/multidomain/embedded/mixeddimensionglue.hh index b1f5f07bcace723bb36e6b647b2d7039514b1f16..231a4bd981b6b6fcfe2692be81014b6c74fa3cb3 100644 --- a/dumux/multidomain/embedded/mixeddimensionglue.hh +++ b/dumux/multidomain/embedded/mixeddimensionglue.hh @@ -25,6 +25,8 @@ #ifndef DUMUX_MULTIDOMAIN_MIXEDDIMENSION_GLUE_HH #define DUMUX_MULTIDOMAIN_MIXEDDIMENSION_GLUE_HH +#warning "This header is deprecated and will be removed after release 3.1. Use multidomain/glue.hh" + #include <dumux/multidomain/glue.hh> #include <dune/grid/common/mcmgmapper.hh> @@ -34,44 +36,8 @@ namespace Dumux { template<class BulkGridView, class LowDimGridView, class BulkMapper = Dune::MultipleCodimMultipleGeomTypeMapper<BulkGridView>, class LowDimMapper = Dune::MultipleCodimMultipleGeomTypeMapper<LowDimGridView>> -class [[deprecated("Use MultiDomainGlue instead. Will be removed after 3.1!")]] -MixedDimensionGlue -: public MultiDomainGlue<LowDimGridView, BulkGridView, LowDimMapper, BulkMapper> -{ - using ParentType = MultiDomainGlue<LowDimGridView, BulkGridView, LowDimMapper, BulkMapper>; -public: - //! Default constructor - MixedDimensionGlue() = default; - - /*! - * \brief constructor from bounding box trees - * \note The definition of Glue::Intersection::inside() has changed! - * inside() used to return the lower-dimensional entities - * (i.e. from LowDimGridView). Now, inside() returns the entities - * of the template argument DomainGridView, i.e. the Dune::GridView - * given as the first template argument. In order to be backwards - * compatible, we flip the order here (and in the definition of ParentType). - */ - template<class BulkTree, class LowDimTree> - MixedDimensionGlue(const BulkTree& bulkTree, const LowDimTree& lowDimTree) - : ParentType(lowDimTree, bulkTree) - {} - - /*! - * \brief construction from bounding box trees - * \note The definition of Glue::Intersection::inside() has changed! - * inside() used to return the lower-dimensional entities - * (i.e. from LowDimGridView). Now, inside() returns the entities - * of the template argument DomainGridView, i.e. the Dune::GridView - * given as the first template argument. In order to be backwards - * compatible, we flip the order here (and in the definition of ParentType). - */ - template<class BulkTree, class LowDimTree> - void build(const BulkTree& bulkTree, const LowDimTree& lowDimTree) - { - ParentType::build(lowDimTree, bulkTree); - } -}; +using MixedDimensionGlue [[deprecated("Use MultiDomainGlue instead. Will be removed after 3.1!")]] + = MultiDomainGlue<BulkGridView, LowDimGridView, BulkMapper, LowDimMapper>; } // end namespace Dumux diff --git a/dumux/multidomain/glue.hh b/dumux/multidomain/glue.hh index 47c96a5cbc7d8f4ed080eb747cefd5f85c9dfa1a..5d0d3c089690097bf6c54a46e0a09a7ec6dbc29a 100644 --- a/dumux/multidomain/glue.hh +++ b/dumux/multidomain/glue.hh @@ -30,10 +30,14 @@ #include <fstream> #include <string> #include <utility> +#include <type_traits> +#include <tuple> +#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 <dune/grid/common/mcmgmapper.hh> @@ -91,6 +95,13 @@ class Intersection using Geometry = Dune::AffineGeometry<ctype, dimIs, dimWorld>; using GlobalPosition = Dune::FieldVector<ctype, dimWorld>; + // 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: Intersection(const DomainTree& domainTree, const TargetTree& targetTree) : domainTree_(domainTree) @@ -107,28 +118,58 @@ public: //! add a pair of neighbor elements void addNeighbors(std::size_t domain, std::size_t target) { - neighbors_[0].push_back(domain); - neighbors_[1].push_back(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 + DomainElement domainEntity(unsigned int n = 0) const + { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); } + + //! get the nth target neighbor entity + TargetElement targetEntity(unsigned int n = 0) const + { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); } + //! get the number of neigbors + [[deprecated("neighbor is deprecated and will be removed after release 3.1. Use numDomainNeighbors() / numTargetNeighbors()")]] std::size_t neighbor(unsigned int side) const - { return neighbors_[side].size(); } + { return std::max(std::get<domainIdx>(neighbors_).size(), std::get<targetIdx>(neighbors_).size()); } //! get the inside (domain) neighbor - DomainElement inside(unsigned int n) const - { return domainTree_.entitySet().entity(neighbors_[0][n]); } + [[deprecated("outside is deprecated and will be removed after release 3.1. Use domainIdx(uint)")]] + DomainElement outside(unsigned int n) const + { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); } //! get the outside (target) neighbor - TargetElement outside(unsigned int n) const - { return targetTree_.entitySet().entity(neighbors_[1][n]); } + [[deprecated("inside is deprecated and will be removed after release 3.1. Use targetIdx(uint)")]] + TargetElement inside(unsigned int n) const + { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); } private: - std::array<std::vector<std::size_t>, 2> neighbors_; + IndexStorage neighbors_; std::vector<GlobalPosition> corners_; const DomainTree& domainTree_;