Commit efd15461 authored by Timo Koch's avatar Timo Koch
Browse files

[md][glue] Fixes change that changed the orientation of inside/outside

To make the interface more clear this commit introduces a new
intersection interface and deprecates the old one.

The old inside(n) is replaced by targetEntity(n)
The old outide(n) is replaced by domainEntity(n)
The old neighbor() is replaces by
numTargetNeighbors()/numDomainNeighbors().

Furthermore num{}Neighbors returns 1 if the {}-dimension is
smaller than the other dimension.
parent 9f8b4305
......@@ -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();
......
......@@ -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>
......
......@@ -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
......
......@@ -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_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment