diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index 4a69f5d7e0a7afea4e82d993f82ba21c9074767e..505a44dfc650c5ced65c6060d6864da6ab8440ee 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -483,14 +483,9 @@ public:
std::vector<PointSource> sources;
asImp_().addPointSources(sources);
- // if there are point sources compute the DOF to point source map
+ // if there are point sources calculate point source locations and save them in a map
if (!sources.empty())
- {
- // calculate point source locations and save them in a map
- PointSourceHelper::computePointSourceMap(*gridGeometry_,
- sources,
- pointSourceMap_);
- }
+ PointSourceHelper::computePointSourceMap(*gridGeometry_, sources, pointSourceMap_, paramGroup());
}
/*!
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 1bd426ad66b90f5b9a69cc02a9f5fa49be678fe2..9633fa6af742e700f707ce8d52b2b0b8d821a2bc 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -28,6 +28,7 @@
#include <functional>
+#include <dune/common/reservedvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/geometry/boundingboxtree.hh>
@@ -278,23 +279,24 @@ public:
//! calculate a DOF index to point source map from given vector of point sources
template<class GridGeometry, class PointSource, class PointSourceMap>
static void computePointSourceMap(const GridGeometry& gridGeometry,
- std::vector<PointSource>& sources,
- PointSourceMap& pointSourceMap)
+ const std::vector<PointSource>& sources,
+ PointSourceMap& pointSourceMap,
+ const std::string& paramGroup = "")
{
- constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
-
const auto& boundingBoxTree = gridGeometry.boundingBoxTree();
- for (auto&& source : sources)
+ for (const auto& s : sources)
{
+ // make local copy of point source for the map
+ auto source = s;
// compute in which elements the point source falls
const auto entities = intersectingEntities(source.position(), boundingBoxTree);
// split the source values equally among all concerned entities
source.setEmbeddings(entities.size()*source.embeddings());
// loop over all concernes elements
- for (unsigned int eIdx : entities)
+ for (const auto eIdx : entities)
{
- if(isBox)
+ if constexpr (GridGeometry::discMethod == DiscretizationMethod::box)
{
// check in which subcontrolvolume(s) we are
const auto element = boundingBoxTree.entitySet().entity(eIdx);
@@ -303,21 +305,21 @@ public:
const auto globalPos = source.position();
// loop over all sub control volumes and check if the point source is inside
- std::vector<unsigned int> scvIndices;
- for (auto&& scv : scvs(fvGeometry))
- {
+ constexpr int dim = GridGeometry::GridView::dimension;
+ Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
+ for (const auto& scv : scvs(fvGeometry))
if (intersectsPointGeometry(globalPos, scv.geometry()))
scvIndices.push_back(scv.indexInElement());
- }
- // for all scvs that where tested positiv add the point sources
+
+ // for all scvs that tested positive add the point sources
// to the element/scv to point source map
- for (auto scvIdx : scvIndices)
+ for (const auto scvIdx : scvIndices)
{
const auto key = std::make_pair(eIdx, scvIdx);
if (pointSourceMap.count(key))
- pointSourceMap.at(key).push_back(source);
+ pointSourceMap.at(key).emplace_back(std::move(source));
else
- pointSourceMap.insert({key, {source}});
+ pointSourceMap.insert({key, {std::move(source)}});
// split equally on the number of matched scvs
auto& s = pointSourceMap.at(key).back();
s.setEmbeddings(scvIndices.size()*s.embeddings());
@@ -328,9 +330,9 @@ public:
// add the pointsource to the DOF map
const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
if (pointSourceMap.count(key))
- pointSourceMap.at(key).push_back(source);
+ pointSourceMap.at(key).emplace_back(std::move(source));
else
- pointSourceMap.insert({key, {source}});
+ pointSourceMap.insert({key, {std::move(source)}});
}
}
}
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh
index 4f495a88cd466466c8c25ae59a549554086bbe66..0b9c09038ccae294cb8eee4ce11746cde9868b9f 100644
--- a/dumux/multidomain/embedded/couplingmanager1d3d.hh
+++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh
@@ -32,6 +32,7 @@
#include <dune/geometry/quadraturerules.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/indextraits.hh>
#include <dumux/multidomain/embedded/pointsourcedata.hh>
#include <dumux/multidomain/embedded/integrationpointsource.hh>
#include <dumux/multidomain/embedded/couplingmanagerbase.hh>
@@ -103,6 +104,7 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::line>
template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
public:
static constexpr EmbeddedCouplingMode couplingMode = EmbeddedCouplingMode::line;
@@ -122,6 +124,9 @@ public:
{
// resize the storage vector
lowDimVolumeInBulkElement_.resize(this->gridView(bulkIdx).size(0));
+ // get references to the grid geometries
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
// compute the low dim volume fractions
for (const auto& is : intersections(this->glue()))
@@ -129,14 +134,14 @@ public:
// all inside elements are identical...
const auto& inside = is.targetEntity(0);
const auto intersectionGeometry = is.geometry();
- const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
- for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
+ for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
}
}
@@ -203,6 +208,7 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::average>
template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
@@ -271,9 +277,11 @@ public:
void computePointSourceData(std::size_t order = 1, bool verbose = false)
{
// Initialize the bulk bounding box tree
- const auto& bulkTree = this->problem(bulkIdx).gridGeometry().boundingBoxTree();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkTree = bulkGridGeometry.boundingBoxTree();
- // initilize the maps
+ // initialize the maps
// do some logging and profiling
Dune::Timer watch;
std::cout << "Initializing the point sources..." << std::endl;
@@ -281,11 +289,11 @@ public:
// clear all internal members like pointsource vectors and stencils
// initializes the point source id counter
this->clear();
- extendedSourceStencil_.stencil().clear();
+ extendedSourceStencil_.clear();
// precompute the vertex indices for efficiency
- this->preComputeVertexIndices(bulkIdx);
- this->preComputeVertexIndices(lowDimIdx);
+ this->precomputeVertexIndices(bulkIdx);
+ this->precomputeVertexIndices(lowDimIdx);
// iterate over all lowdim elements
const auto& lowDimProblem = this->problem(lowDimIdx);
@@ -295,7 +303,7 @@ public:
const auto lowDimGeometry = lowDimElement.geometry();
const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(lowDimGeometry.type(), order);
- const auto lowDimElementIdx = lowDimProblem.gridGeometry().elementMapper().index(lowDimElement);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement);
// apply the Gaussian quadrature rule and define point sources at each quadrature point
// note that the approximation is not optimal if
@@ -323,36 +331,43 @@ public:
static const auto numIp = getParam<int>("MixedDimension.NumCircleSegments");
const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx);
const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0);
- const auto weight = 2*M_PI*radius/numIp;
+ const auto circleAvgWeight = 2*M_PI*radius/numIp;
const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp);
- std::vector<Scalar> circleIpWeight(circlePoints.size());
- std::vector<std::size_t> circleStencil(circlePoints.size());
+ std::vector<Scalar> circleIpWeight; circleIpWeight.reserve(circlePoints.size());
+ std::vector<GridIndex<bulkIdx>> circleStencil; circleStencil.reserve(circlePoints.size());
+
// for box
- std::unordered_map<std::size_t, std::vector<std::size_t> > circleCornerIndices;
+ std::vector<const std::vector<GridIndex<bulkIdx>>*> circleCornerIndices;
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
- std::unordered_map<std::size_t, ShapeValues> circleShapeValues;
+ std::vector<ShapeValues> circleShapeValues;
+ // go over all points of the average operator
+ int insideCirclePoints = 0;
for (int k = 0; k < circlePoints.size(); ++k)
{
const auto circleBulkElementIndices = intersectingEntities(circlePoints[k], bulkTree);
if (circleBulkElementIndices.empty())
continue;
- const auto bulkElementIdx = circleBulkElementIndices[0];
- circleStencil[k] = bulkElementIdx;
- circleIpWeight[k] = weight;
-
- if (isBox<bulkIdx>())
+ ++insideCirclePoints;
+ const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices.size();
+ for (const auto bulkElementIdx : circleBulkElementIndices)
{
- if (!static_cast<bool>(circleCornerIndices.count(bulkElementIdx)))
+ circleStencil.push_back(bulkElementIdx);
+ circleIpWeight.push_back(localCircleAvgWeight);
+
+ // precompute interpolation data for box scheme for each cut element
+ if (isBox<bulkIdx>())
{
- const auto bulkElement = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx);
- circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx);
+ const auto bulkElement = bulkGridGeometry.element(bulkElementIdx);
+ circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx)));
// evaluate shape functions at the integration point
const auto bulkGeometry = bulkElement.geometry();
- this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]);
+ ShapeValues shapeValues;
+ this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, circlePoints[k], shapeValues);
+ circleShapeValues.emplace_back(std::move(shapeValues));
}
}
}
@@ -364,7 +379,7 @@ public:
for (const auto& vertices : circleCornerIndices)
{
this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(),
- vertices.second.begin(), vertices.second.end());
+ vertices->begin(), vertices->end());
}
}
@@ -390,10 +405,10 @@ public:
// the actual solution dependent source term
PointSourceData psData;
- if (isBox<lowDimIdx>())
+ if constexpr (isBox<lowDimIdx>())
{
ShapeValues shapeValues;
- this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
+ this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues);
psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
}
else
@@ -402,13 +417,13 @@ public:
}
// add data needed to compute integral over the circle
- if (isBox<bulkIdx>())
+ if constexpr (isBox<bulkIdx>())
{
psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil);
- const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
+ const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry();
ShapeValues shapeValues;
- this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues);
+ this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues);
psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
}
else
@@ -440,7 +455,7 @@ public:
for (const auto& vertices : circleCornerIndices)
{
extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(),
- vertices.second.begin(), vertices.second.end());
+ vertices->begin(), vertices->end());
}
}
@@ -495,6 +510,9 @@ public:
{
// resize the storage vector
lowDimVolumeInBulkElement_.resize(this->gridView(bulkIdx).size(0));
+ // get references to the grid geometries
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
// compute the low dim volume fractions
for (const auto& is : intersections(this->glue()))
@@ -502,14 +520,14 @@ public:
// all inside elements are identical...
const auto& inside = is.targetEntity(0);
const auto intersectionGeometry = is.geometry();
- const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
- for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
+ for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
}
}
@@ -580,6 +598,7 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::cylindersource
template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
@@ -617,9 +636,11 @@ public:
void computePointSourceData(std::size_t order = 1, bool verbose = false)
{
// Initialize the bulk bounding box tree
- const auto& bulkTree = this->problem(bulkIdx).gridGeometry().boundingBoxTree();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkTree = bulkGridGeometry.boundingBoxTree();
- // initilize the maps
+ // initialize the maps
// do some logging and profiling
Dune::Timer watch;
std::cout << "Initializing the point sources..." << std::endl;
@@ -629,8 +650,8 @@ public:
this->clear();
// precompute the vertex indices for efficiency
- this->preComputeVertexIndices(bulkIdx);
- this->preComputeVertexIndices(lowDimIdx);
+ this->precomputeVertexIndices(bulkIdx);
+ this->precomputeVertexIndices(lowDimIdx);
// iterate over all lowdim elements
const auto& lowDimProblem = this->problem(lowDimIdx);
@@ -640,7 +661,7 @@ public:
const auto lowDimGeometry = lowDimElement.geometry();
const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(lowDimGeometry.type(), order);
- const auto lowDimElementIdx = lowDimProblem.gridGeometry().elementMapper().index(lowDimElement);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement);
// apply the Gaussian quadrature rule and define point sources at each quadrature point
// note that the approximation is not optimal if
@@ -693,11 +714,11 @@ public:
// the actual solution dependent source term
PointSourceData psData;
- if (isBox<lowDimIdx>())
+ if constexpr (isBox<lowDimIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
ShapeValues shapeValues;
- this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
+ this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues);
psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
}
else
@@ -706,12 +727,12 @@ public:
}
// add data needed to compute integral over the circle
- if (isBox<bulkIdx>())
+ if constexpr (isBox<bulkIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
- const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
+ const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry();
ShapeValues shapeValues;
- this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, circlePos, shapeValues);
+ this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, circlePos, shapeValues);
psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
}
else
@@ -773,6 +794,9 @@ public:
{
// resize the storage vector
lowDimVolumeInBulkElement_.resize(this->gridView(bulkIdx).size(0));
+ // get references to the grid geometries
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
// compute the low dim volume fractions
for (const auto& is : intersections(this->glue()))
@@ -780,14 +804,14 @@ public:
// all inside elements are identical...
const auto& inside = is.targetEntity(0);
const auto intersectionGeometry = is.geometry();
- const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
- for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
+ for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
}
}
@@ -856,6 +880,7 @@ class EmbeddedCouplingManager1d3d<MDTraits, EmbeddedCouplingMode::kernel>
template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
@@ -921,7 +946,7 @@ public:
*/
void computePointSourceData(std::size_t order = 1, bool verbose = false)
{
- // initilize the maps
+ // initialize the maps
// do some logging and profiling
Dune::Timer watch;
std::cout << "Initializing the point sources..." << std::endl;
@@ -931,14 +956,14 @@ public:
this->clear();
bulkSourceIds_.clear();
bulkSourceWeights_.clear();
- extendedSourceStencil_.stencil().clear();
+ extendedSourceStencil_.clear();
// precompute the vertex indices for efficiency for the box method
- this->preComputeVertexIndices(bulkIdx);
- this->preComputeVertexIndices(lowDimIdx);
+ this->precomputeVertexIndices(bulkIdx);
+ this->precomputeVertexIndices(lowDimIdx);
- const auto& bulkFvGridGeometry = this->problem(bulkIdx).gridGeometry();
- const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
bulkSourceIds_.resize(this->gridView(bulkIdx).size(0));
bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0));
@@ -958,7 +983,7 @@ public:
// get the Gaussian quadrature rule for the local intersection
const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
- const std::size_t lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(inside);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// iterate over all quadrature points and place a source
// for 1d: make a new point source
@@ -966,10 +991,10 @@ public:
for (auto&& qp : quad)
{
// compute the coupling stencils
- for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
+ for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside);
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
// each quadrature point will be a point source for the sub problem
const auto globalPos = intersectionGeometry.global(qp.position());
@@ -984,12 +1009,12 @@ public:
// the actual solution dependent source term
PointSourceData psData;
- if (isBox<lowDimIdx>())
+ if constexpr (isBox<lowDimIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
- const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
+ const auto lowDimGeometry = lowDimGridGeometry.element(lowDimElementIdx).geometry();
ShapeValues shapeValues;
- this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
+ this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues);
psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
}
else
@@ -998,12 +1023,12 @@ public:
}
// add data needed to compute integral over the circle
- if (isBox<bulkIdx>())
+ if constexpr (isBox<bulkIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
- const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
+ const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry();
ShapeValues shapeValues;
- this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues);
+ this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues);
psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
}
else
@@ -1078,6 +1103,9 @@ public:
{
// resize the storage vector
lowDimVolumeInBulkElement_.resize(this->gridView(bulkIdx).size(0));
+ // get references to the grid geometries
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
// compute the low dim volume fractions
for (const auto& is : intersections(this->glue()))
@@ -1085,14 +1113,14 @@ public:
// all inside elements are identical...
const auto& inside = is.targetEntity(0);
const auto intersectionGeometry = is.geometry();
- const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
+ const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
- for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
+ for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
}
}
@@ -1127,11 +1155,11 @@ public:
}
//! return all source ids for a bulk elements
- const std::vector<std::size_t> bulkSourceIds(std::size_t eIdx) const
+ const std::vector<std::size_t> bulkSourceIds(GridIndex<bulkIdx> eIdx) const
{ return bulkSourceIds_[eIdx]; }
//! return all source ids for a bulk elements
- const std::vector<Scalar> bulkSourceWeights(std::size_t eIdx) const
+ const std::vector<Scalar> bulkSourceWeights(GridIndex<bulkIdx> eIdx) const
{ return bulkSourceWeights_[eIdx]; }
// \}
@@ -1139,7 +1167,7 @@ public:
private:
//! TODO: How to optimize this?
void computeBulkSource(const GlobalPosition& globalPos, const Scalar kernelWidth,
- std::size_t id, std::size_t lowDimElementIdx, std::size_t coupledBulkElementIdx,
+ std::size_t id, GridIndex<lowDimIdx> lowDimElementIdx, GridIndex<bulkIdx> coupledBulkElementIdx,
Scalar pointSourceWeight)
{
// make sure it is mass conservative
@@ -1162,7 +1190,8 @@ private:
if (weight > 1e-13)
{
- const auto bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(element);
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(element);
bulkSourceIds_[bulkElementIdx].push_back(id);
bulkSourceWeights_[bulkElementIdx].push_back(weight*pointSourceWeight);
mask[bulkElementIdx] = true;
@@ -1188,7 +1217,7 @@ private:
}
}
- for (std::size_t eIdx = 0; eIdx < bulkSourceWeights_.size(); ++eIdx)
+ for (GridIndex<bulkIdx> eIdx = 0; eIdx < bulkSourceWeights_.size(); ++eIdx)
if (mask[eIdx]) bulkSourceWeights_[eIdx].back() /= checkSum;
// balance error of the quadrature rule -> TODO: what to do at boundaries
diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh
index 678eb0ad4ce9b69dbd97a6d3f1475a0564f1f15b..b130b9edd1021e8fbe7d2675fa4e68c3b90dbafb 100644
--- a/dumux/multidomain/embedded/couplingmanagerbase.hh
+++ b/dumux/multidomain/embedded/couplingmanagerbase.hh
@@ -91,18 +91,17 @@ class EmbeddedCouplingManagerBase
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using ElementMapper = typename GridGeometry<id>::ElementMapper;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
+ template<std::size_t id> using CouplingStencil = std::vector<GridIndex<id>>;
- enum {
- bulkDim = GridView<bulkIdx>::dimension,
- lowDimDim = GridView<lowDimIdx>::dimension,
- dimWorld = GridView<bulkIdx>::dimensionworld
- };
+ static constexpr int bulkDim = GridView<bulkIdx>::dimension;
+ static constexpr int lowDimDim = GridView<lowDimIdx>::dimension;
+ static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld;
template<std::size_t id>
static constexpr bool isBox()
{ return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
- using CouplingStencil = std::vector<std::size_t>;
using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
using GlueType = MultiDomainGlue<GridView<bulkIdx>, GridView<lowDimIdx>, ElementMapper<bulkIdx>, ElementMapper<lowDimIdx>>;
@@ -112,13 +111,13 @@ public:
//! export the point source traits
using PointSourceTraits = PSTraits;
//! export stencil types
- using CouplingStencils = std::unordered_map<std::size_t, CouplingStencil>;
+ template<std::size_t id> using CouplingStencils = std::unordered_map<GridIndex<id>, CouplingStencil<id>>;
/*!
* \brief call this after grid adaption
*/
- void updateAfterGridAdaption(std::shared_ptr<const GridGeometry<bulkIdx>> bulkFvGridGeometry,
- std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimFvGridGeometry)
+ void updateAfterGridAdaption(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry,
+ std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry)
{
glue_ = std::make_shared<GlueType>();
}
@@ -126,10 +125,10 @@ public:
/*!
* \brief Constructor
*/
- EmbeddedCouplingManagerBase(std::shared_ptr<const GridGeometry<bulkIdx>> bulkFvGridGeometry,
- std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimFvGridGeometry)
+ EmbeddedCouplingManagerBase(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry,
+ std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry)
{
- updateAfterGridAdaption(bulkFvGridGeometry, lowDimFvGridGeometry);
+ updateAfterGridAdaption(bulkGridGeometry, lowDimGridGeometry);
}
/*!
@@ -170,9 +169,9 @@ public:
* \note This function has to be implemented by all coupling managers for all combinations of i and j
*/
template<std::size_t i, std::size_t j>
- const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
- const Element<i>& element,
- Dune::index_constant<j> domainJ) const
+ const CouplingStencil<j>& couplingStencil(Dune::index_constant<i> domainI,
+ const Element<i>& element,
+ Dune::index_constant<j> domainJ) const
{
static_assert(i != j, "A domain cannot be coupled to itself!");
@@ -180,7 +179,7 @@ public:
if (couplingStencils(domainI).count(eIdx))
return couplingStencils(domainI).at(eIdx);
else
- return emptyStencil_;
+ return emptyStencil(domainI);
}
/*!
@@ -243,11 +242,11 @@ public:
clear();
// precompute the vertex indices for efficiency for the box method
- this->preComputeVertexIndices(bulkIdx);
- this->preComputeVertexIndices(lowDimIdx);
+ this->precomputeVertexIndices(bulkIdx);
+ this->precomputeVertexIndices(lowDimIdx);
- const auto& bulkFvGridGeometry = this->problem(bulkIdx).gridGeometry();
- const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
// intersect the bounding box trees
glueGrids();
@@ -263,7 +262,7 @@ public:
// get the Gaussian quadrature rule for the local intersection
const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
- const std::size_t lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(inside);
+ const std::size_t lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
// iterate over all quadrature points
for (auto&& qp : quad)
@@ -272,7 +271,7 @@ public:
for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
{
const auto& outside = is.domainEntity(outsideIdx);
- const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside);
+ const std::size_t bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
// each quadrature point will be a point source for the sub problem
const auto globalPos = intersectionGeometry.global(qp.position());
@@ -288,7 +287,7 @@ public:
// the actual solution dependent source term
PointSourceData psData;
- if (isBox<lowDimIdx>())
+ if constexpr (isBox<lowDimIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
@@ -302,7 +301,7 @@ public:
}
// add data needed to compute integral over the circle
- if (isBox<bulkIdx>())
+ if constexpr (isBox<bulkIdx>())
{
using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
@@ -381,17 +380,11 @@ public:
//! Return data for a bulk point source with the identifier id
PrimaryVariables<bulkIdx> bulkPriVars(std::size_t id) const
- {
- auto& data = pointSourceData_[id];
- return data.interpolateBulk(this->curSol()[bulkIdx]);
- }
+ { return pointSourceData_[id].interpolateBulk(this->curSol()[bulkIdx]); }
//! Return data for a low dim point source with the identifier id
PrimaryVariables<lowDimIdx> lowDimPriVars(std::size_t id) const
- {
- auto& data = pointSourceData_[id];
- return data.interpolateLowDim(this->curSol()[lowDimIdx]);
- }
+ { return pointSourceData_[id].interpolateLowDim(this->curSol()[lowDimIdx]); }
//! return the average distance to the coupled bulk cell center
Scalar averageDistance(std::size_t id) const
@@ -412,25 +405,26 @@ public:
//! Return reference to bulk coupling stencil member of domain i
template<std::size_t i>
- const CouplingStencils& couplingStencils(Dune::index_constant<i> dom) const
- { return (i == bulkIdx) ? bulkCouplingStencils_ : lowDimCouplingStencils_; }
+ const CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom) const
+ { return std::get<i>(couplingStencils_); }
//! Return reference to point source data vector member
const std::vector<PointSourceData>& pointSourceData() const
{ return pointSourceData_; }
//! Return a reference to an empty stencil
- const CouplingStencil& emptyStencil() const
- { return emptyStencil_; }
+ template<std::size_t i>
+ const CouplingStencil<i>& emptyStencil(Dune::index_constant<i> dom) const
+ { return std::get<i>(emptyStencil_); }
protected:
//! computes the vertex indices per element for the box method
template<std::size_t id>
- void preComputeVertexIndices(Dune::index_constant<id> domainIdx)
+ void precomputeVertexIndices(Dune::index_constant<id> domainIdx)
{
// fill helper structure for box discretization
- if (isBox<domainIdx>())
+ if constexpr (isBox<domainIdx>())
{
this->vertexIndices(domainIdx).resize(gridView(domainIdx).size(0));
for (const auto& element : elements(gridView(domainIdx)))
@@ -445,19 +439,17 @@ protected:
}
//! compute the shape function for a given point and geometry
- template<std::size_t i, class FVGG, class Geometry, class ShapeValues, typename std::enable_if_t<FVGG::discMethod == DiscretizationMethod::box, int> = 0>
+ template<std::size_t i, class FVGG, class Geometry, class ShapeValues>
void getShapeValues(Dune::index_constant<i> domainI, const FVGG& gridGeometry, const Geometry& geo, const GlobalPosition& globalPos, ShapeValues& shapeValues)
{
- const auto ipLocal = geo.local(globalPos);
- const auto& localBasis = this->problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis();
- localBasis.evaluateFunction(ipLocal, shapeValues);
- }
-
- //! compute the shape function for a given point and geometry
- template<std::size_t i, class FVGG, class Geometry, class ShapeValues, typename std::enable_if_t<FVGG::discMethod != DiscretizationMethod::box, int> = 0>
- void getShapeValues(Dune::index_constant<i> domainI, const FVGG& gridGeometry, const Geometry& geo, const GlobalPosition& globalPos, ShapeValues& shapeValues)
- {
- DUNE_THROW(Dune::InvalidStateException, "Shape values requested for other discretization than box!");
+ if constexpr (FVGG::discMethod == DiscretizationMethod::box)
+ {
+ const auto ipLocal = geo.local(globalPos);
+ const auto& localBasis = this->problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis();
+ localBasis.evaluateFunction(ipLocal, shapeValues);
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Shape values requested for other discretization than box!");
}
//! Clear all internal data members
@@ -465,11 +457,11 @@ protected:
{
pointSources(bulkIdx).clear();
pointSources(lowDimIdx).clear();
+ couplingStencils(bulkIdx).clear();
+ couplingStencils(lowDimIdx).clear();
+ vertexIndices(bulkIdx).clear();
+ vertexIndices(lowDimIdx).clear();
pointSourceData_.clear();
- bulkCouplingStencils_.clear();
- lowDimCouplingStencils_.clear();
- bulkVertexIndices_.clear();
- lowDimVertexIndices_.clear();
averageDistanceToBulkCell_.clear();
idCounter_ = 0;
@@ -478,23 +470,20 @@ protected:
//! compute the intersections between the two grids
void glueGrids()
{
- const auto& bulkFvGridGeometry = this->problem(bulkIdx).gridGeometry();
- const auto& lowDimFvGridGeometry = this->problem(lowDimIdx).gridGeometry();
+ const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
+ const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
// intersect the bounding box trees
- glue_->build(bulkFvGridGeometry.boundingBoxTree(), lowDimFvGridGeometry.boundingBoxTree());
+ glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree());
}
template<class Geometry, class GlobalPosition>
- Scalar computeDistance(const Geometry& geometry, const GlobalPosition& p)
+ Scalar computeDistance(const Geometry& geometry, const GlobalPosition& p) const
{
Scalar avgDist = 0.0;
const auto& quad = Dune::QuadratureRules<Scalar, bulkDim>::rule(geometry.type(), 5);
for (auto&& qp : quad)
- {
- const auto globalPos = geometry.global(qp.position());
- avgDist += (globalPos-p).two_norm()*qp.weight();
- }
+ avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight();
return avgDist;
}
@@ -513,18 +502,18 @@ protected:
//! Return reference to bulk coupling stencil member of domain i
template<std::size_t i>
- CouplingStencils& couplingStencils(Dune::index_constant<i> dom)
- { return (i == 0) ? bulkCouplingStencils_ : lowDimCouplingStencils_; }
+ CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom)
+ { return std::get<i>(couplingStencils_); }
//! Return a reference to the vertex indices
template<std::size_t i>
- std::vector<std::size_t>& vertexIndices(Dune::index_constant<i> dom, std::size_t eIdx)
- { return (i == 0) ? bulkVertexIndices_[eIdx] : lowDimVertexIndices_[eIdx]; }
+ std::vector<GridIndex<i>>& vertexIndices(Dune::index_constant<i> dom, GridIndex<i> eIdx)
+ { return std::get<i>(vertexIndices_)[eIdx]; }
//! Return a reference to the vertex indices container
template<std::size_t i>
- std::vector<std::vector<std::size_t>>& vertexIndices(Dune::index_constant<i> dom)
- { return (i == 0) ? bulkVertexIndices_ : lowDimVertexIndices_; }
+ std::vector<std::vector<GridIndex<i>>>& vertexIndices(Dune::index_constant<i> dom)
+ { return std::get<i>(vertexIndices_); }
const GlueType& glue() const
{ return *glue_; }
@@ -544,15 +533,14 @@ private:
//! the point source in both domains
std::tuple<std::vector<PointSource<bulkIdx>>, std::vector<PointSource<lowDimIdx>>> pointSources_;
- mutable std::vector<PointSourceData> pointSourceData_;
+ std::vector<PointSourceData> pointSourceData_;
std::vector<Scalar> averageDistanceToBulkCell_;
//! Stencil data
- std::vector<std::vector<std::size_t>> bulkVertexIndices_;
- std::vector<std::vector<std::size_t>> lowDimVertexIndices_;
- CouplingStencils bulkCouplingStencils_;
- CouplingStencils lowDimCouplingStencils_;
- CouplingStencil emptyStencil_;
+ std::tuple<std::vector<std::vector<GridIndex<bulkIdx>>>,
+ std::vector<std::vector<GridIndex<lowDimIdx>>>> vertexIndices_;
+ std::tuple<CouplingStencils<bulkIdx>, CouplingStencils<lowDimIdx>> couplingStencils_;
+ std::tuple<CouplingStencil<bulkIdx>, CouplingStencil<lowDimIdx>> emptyStencil_;
//! The glue object
std::shared_ptr<GlueType> glue_;
diff --git a/dumux/multidomain/embedded/extendedsourcestencil.hh b/dumux/multidomain/embedded/extendedsourcestencil.hh
index 9d1baae3e6ae44e23a343dfa0f0de1f28ec924db..1e1a3afa62a244742dd4773a430d8040ad6ff87e 100644
--- a/dumux/multidomain/embedded/extendedsourcestencil.hh
+++ b/dumux/multidomain/embedded/extendedsourcestencil.hh
@@ -34,8 +34,7 @@
#include <dumux/common/numericdifferentiation.hh>
#include <dumux/discretization/method.hh>
-namespace Dumux {
-namespace EmbeddedCoupling {
+namespace Dumux::EmbeddedCoupling {
/*!
* \ingroup EmbeddedCoupling
@@ -53,8 +52,10 @@ class ExtendedSourceStencil
template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
- static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
- static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
+ template<std::size_t id>
+ static constexpr auto subDomainIdx = typename MDTraits::template SubDomain<id>::Index();
+ static constexpr auto bulkIdx = subDomainIdx<0>;
+ static constexpr auto lowDimIdx = subDomainIdx<1>;
template<std::size_t id>
static constexpr bool isBox()
@@ -73,8 +74,7 @@ public:
for (const auto& element : elements(couplingManager.gridView(domainI)))
{
const auto& dofs = extendedSourceStencil_(couplingManager, domainI, element);
-
- if (isBox<domainI>())
+ if constexpr (isBox<domainI>())
{
for (int i = 0; i < element.subEntities(GridView<domainI>::dimension); ++i)
for (const auto globalJ : dofs)
@@ -104,7 +104,7 @@ public:
JacobianMatrixDiagBlock& A,
GridVariables& gridVariables) const
{
- constexpr auto numEq = std::decay_t<decltype(curSol[domainI][0])>::dimension;
+ constexpr auto numEq = std::decay_t<decltype(curSol[domainI][0])>::size();
const auto& elementI = localAssemblerI.element();
// only do something if we have an extended stencil
@@ -141,7 +141,7 @@ public:
partialDerivs, origResidual, numDiffMethod);
// update the global stiffness matrix with the current partial derivatives
- for (auto&& scvJ : scvs(localAssemblerI.fvGeometry()))
+ for (const auto& scvJ : scvs(localAssemblerI.fvGeometry()))
{
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
{
@@ -159,30 +159,32 @@ public:
}
}
+ //! clear the internal data
+ void clear() { sourceStencils_.clear(); }
+
//! return a reference to the stencil
- typename CouplingManager::CouplingStencils& stencil()
+ typename CouplingManager::template CouplingStencils<bulkIdx>& stencil()
{ return sourceStencils_; }
private:
- //! the extended source stencil for the bulk domain due to the source average
- const std::vector<std::size_t>& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<0> bulkDomain, const Element<0>& bulkElement) const
+ //! the extended source stencil due to the source average (always empty for lowdim, but may be filled for bulk)
+ template<std::size_t id>
+ const auto& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<id> bulkDomain, const Element<id>& bulkElement) const
{
- const auto bulkElementIdx = couplingManager.problem(bulkIdx).gridGeometry().elementMapper().index(bulkElement);
- if (sourceStencils_.count(bulkElementIdx))
- return sourceStencils_.at(bulkElementIdx);
- else
- return couplingManager.emptyStencil();
- }
+ if constexpr (subDomainIdx<id> == bulkIdx)
+ {
+ const auto bulkElementIdx = couplingManager.problem(bulkIdx).gridGeometry().elementMapper().index(bulkElement);
+ if (sourceStencils_.count(bulkElementIdx))
+ return sourceStencils_.at(bulkElementIdx);
+ }
- //! the extended source stencil for the low dim domain is empty
- const std::vector<std::size_t>& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<1> bulkDomain, const Element<1>& lowDimElement) const
- { return couplingManager.emptyStencil(); }
+ return couplingManager.emptyStencil(subDomainIdx<id>);
+ }
//! the additional stencil for the kernel evaluations / circle averages
- typename CouplingManager::CouplingStencils sourceStencils_;
+ typename CouplingManager::template CouplingStencils<bulkIdx> sourceStencils_;
};
-} // end namespace EmbeddedCoupling
-} // end namespace Dumux
+} // end namespace Dumux::EmbeddedCoupling
#endif
diff --git a/dumux/multidomain/embedded/integrationpointsource.hh b/dumux/multidomain/embedded/integrationpointsource.hh
index 74ac494c9f378c6394ee78000e1155c8c14c6c60..e842481517dd8229578d6cb60fcc9b8809acdb2a 100644
--- a/dumux/multidomain/embedded/integrationpointsource.hh
+++ b/dumux/multidomain/embedded/integrationpointsource.hh
@@ -30,6 +30,7 @@
#include <type_traits>
#include <dune/common/reservedvector.hh>
#include <dumux/common/pointsource.hh>
+#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
namespace Dumux {
@@ -47,37 +48,64 @@ class IntegrationPointSource : public IdPointSource<GlobalPosition, SourceValues
public:
//! Constructor for integration point sources
+ [[deprecated("Call constructor with a single element index instead! Will be removed after 3.3")]]
IntegrationPointSource(GlobalPosition pos, SourceValues values, IdType id,
Scalar qpweight, Scalar integrationElement,
const std::vector<std::size_t>& elementIndices)
: ParentType(pos, values, id),
qpweight_(qpweight), integrationElement_(integrationElement),
- elementIndices_(elementIndices) {}
+ elementIndex_(elementIndices[0]) {}
//! Constructor for integration point sources, when there is no
// value known at the time of initialization
+ [[deprecated("Call constructor with a single element index instead! Will be removed after 3.3")]]
IntegrationPointSource(GlobalPosition pos, IdType id,
Scalar qpweight, Scalar integrationElement,
const std::vector<std::size_t>& elementIndices)
: ParentType(pos, id),
qpweight_(qpweight), integrationElement_(integrationElement),
- elementIndices_(elementIndices) {}
+ elementIndex_(elementIndices[0]) {}
+ //! Constructor for integration point sources
+ IntegrationPointSource(GlobalPosition pos, SourceValues values, IdType id,
+ Scalar qpweight, Scalar integrationElement,
+ std::size_t elementIndex)
+ : ParentType(pos, values, id)
+ , qpweight_(qpweight)
+ , integrationElement_(integrationElement)
+ , elementIndex_(elementIndex)
+ {}
+
+ //! Constructor for integration point sources, when there is no
+ //! value known at the time of initialization
+ IntegrationPointSource(GlobalPosition pos, IdType id,
+ Scalar qpweight, Scalar integrationElement,
+ std::size_t elementIndex)
+ : ParentType(pos, id)
+ , qpweight_(qpweight)
+ , integrationElement_(integrationElement)
+ , elementIndex_(elementIndex)
+ {}
Scalar quadratureWeight() const
- {
- return qpweight_;
- }
+ { return qpweight_; }
Scalar integrationElement() const
- {
- return integrationElement_;
- }
+ { return integrationElement_; }
- const std::vector<std::size_t>& elementIndices() const
- {
- return elementIndices_;
- }
+ void setQuadratureWeight(const Scalar qpWeight)
+ { return qpweight_ = qpWeight; }
+
+ void setIntegrationElement(const Scalar ie)
+ { integrationElement_ = ie; }
+
+ [[deprecated("Use elementIndex() instead. Will be removed after 3.3")]]
+ std::vector<std::size_t> elementIndices() const
+ { return std::vector<std::size_t>({elementIndex_}); }
+
+ //! The index of the element this intergration point source is associated with
+ std::size_t elementIndex() const
+ { return elementIndex_; }
//! Convenience = operator overload modifying only the values
IntegrationPointSource& operator= (const SourceValues& values)
@@ -96,7 +124,7 @@ public:
private:
Scalar qpweight_;
Scalar integrationElement_;
- std::vector<std::size_t> elementIndices_;
+ std::size_t elementIndex_;
};
/*!
@@ -110,29 +138,29 @@ public:
//! calculate a DOF index to point source map from given vector of point sources
template<class GridGeometry, class PointSource, class PointSourceMap>
static void computePointSourceMap(const GridGeometry& gridGeometry,
- std::vector<PointSource>& sources,
- PointSourceMap& pointSourceMap)
+ const std::vector<PointSource>& sources,
+ PointSourceMap& pointSourceMap,
+ const std::string& paramGroup = "")
{
- for (auto&& source : sources)
+ for (const auto& source : sources)
{
- // compute in which elements the point source falls
- const auto& entities = source.elementIndices();
- // split the source values equally among all concerned entities
- source.setEmbeddings(source.embeddings()*entities.size());
- // loop over all intersected elements
- for (unsigned int eIdx : entities)
+ // get the index of the element in which the point source falls
+ const auto eIdx = source.elementIndex();
+ if constexpr (GridGeometry::discMethod == DiscretizationMethod::box)
{
- if (GridGeometry::discMethod == DiscretizationMethod::box)
+ // check in which subcontrolvolume(s) we are
+ const auto element = gridGeometry.boundingBoxTree().entitySet().entity(eIdx);
+ auto fvGeometry = localView(gridGeometry);
+ fvGeometry.bindElement(element);
+ const auto globalPos = source.position();
+
+ static const bool boxPointSourceLumping = getParamFromGroup<bool>(paramGroup, "PointSource.EnableBoxLumping", true);
+ if (boxPointSourceLumping)
{
- // check in which subcontrolvolume(s) we are
- const auto element = gridGeometry.boundingBoxTree().entitySet().entity(eIdx);
- auto fvGeometry = localView(gridGeometry);
- fvGeometry.bindElement(element);
- const auto globalPos = source.position();
// loop over all sub control volumes and check if the point source is inside
constexpr int dim = GridGeometry::GridView::dimension;
Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
- for (auto&& scv : scvs(fvGeometry))
+ for (const auto& scv : scvs(fvGeometry))
if (intersectsPointGeometry(globalPos, scv.geometry()))
scvIndices.push_back(scv.indexInElement());
@@ -150,16 +178,38 @@ public:
s.setEmbeddings(scvIndices.size()*s.embeddings());
}
}
+
+ // distribute the sources using the basis function weights
else
{
- // add the pointsource to the DOF map
- const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
- if (pointSourceMap.count(key))
- pointSourceMap.at(key).push_back(source);
- else
- pointSourceMap.insert({key, {source}});
+ const auto& localBasis = fvGeometry.feLocalBasis();
+ const auto ipLocal = element.geometry().local(globalPos);
+ using Scalar = std::decay_t<decltype(source.values()[0])>;
+ std::vector<typename Dune::FieldVector<Scalar, 1>> shapeValues;
+ localBasis.evaluateFunction(ipLocal, shapeValues);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto key = std::make_pair(eIdx, scv.indexInElement());
+ if (pointSourceMap.count(key))
+ pointSourceMap.at(key).push_back(source);
+ else
+ pointSourceMap.insert({key, {source}});
+
+ // adjust the integration element
+ auto& s = pointSourceMap.at(key).back();
+ s.setIntegrationElement(shapeValues[scv.indexInElement()]*s.integrationElement());
+ }
}
}
+ else
+ {
+ // add the pointsource to the DOF map
+ const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
+ if (pointSourceMap.count(key))
+ pointSourceMap.at(key).push_back(source);
+ else
+ pointSourceMap.insert({key, {source}});
+ }
}
}
};
diff --git a/dumux/multidomain/embedded/pointsourcedata.hh b/dumux/multidomain/embedded/pointsourcedata.hh
index ef7689e7cf7bdf109b5f139067becc0f63b7674e..d6cbab29645e1b3d5ef5b44645309c76a9820c3f 100644
--- a/dumux/multidomain/embedded/pointsourcedata.hh
+++ b/dumux/multidomain/embedded/pointsourcedata.hh
@@ -22,13 +22,13 @@
* \brief Data associated with a point source
*/
-#ifndef DUMUX_POINTSOURCEDATA_HH
-#define DUMUX_POINTSOURCEDATA_HH
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
#include <vector>
-#include <unordered_map>
#include <dune/common/fvector.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/indextraits.hh>
#include <dumux/discretization/method.hh>
namespace Dumux {
@@ -44,59 +44,60 @@ class PointSourceData
using Scalar = typename MDTraits::Scalar;
using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
- // obtain the type tags of the sub problems
- using BulkTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
- using LowDimTypeTag = typename MDTraits::template SubDomain<1>::TypeTag;
+ template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
+ template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
+ template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
- using BulkPrimaryVariables = GetPropType<BulkTypeTag, Properties::PrimaryVariables>;
- using LowDimPrimaryVariables = GetPropType<LowDimTypeTag, Properties::PrimaryVariables>;
+ static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+ static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
- using BulkSolutionVector = GetPropType<BulkTypeTag, Properties::SolutionVector>;
- using LowDimSolutionVector = GetPropType<LowDimTypeTag, Properties::SolutionVector>;
-
- enum {
- bulkIsBox = GetPropType<BulkTypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box,
- lowDimIsBox = GetPropType<LowDimTypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box
- };
+ template<std::size_t id>
+ static constexpr bool isBox()
+ { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
public:
void addBulkInterpolation(const ShapeValues& shapeValues,
- const std::vector<std::size_t>& cornerIndices,
- std::size_t eIdx)
+ const std::vector<GridIndex<bulkIdx>>& cornerIndices,
+ GridIndex<bulkIdx> eIdx)
{
+ static_assert(isBox<bulkIdx>(), "This interface is only available for the box method.");
bulkElementIdx_ = eIdx;
bulkCornerIndices_ = cornerIndices;
bulkShapeValues_ = shapeValues;
}
void addLowDimInterpolation(const ShapeValues& shapeValues,
- const std::vector<std::size_t>& cornerIndices,
- std::size_t eIdx)
+ const std::vector<GridIndex<lowDimIdx>>& cornerIndices,
+ GridIndex<lowDimIdx> eIdx)
{
+ static_assert(isBox<lowDimIdx>(), "This interface is only available for the box method.");
lowDimElementIdx_ = eIdx;
lowDimCornerIndices_ = cornerIndices;
lowDimShapeValues_ = shapeValues;
}
- void addBulkInterpolation(std::size_t eIdx)
+ void addBulkInterpolation(GridIndex<bulkIdx> eIdx)
{
- assert(!bulkIsBox);
+ static_assert(!isBox<bulkIdx>(), "This interface is not available for the box method.");
bulkElementIdx_ = eIdx;
}
- void addLowDimInterpolation(std::size_t eIdx)
+ void addLowDimInterpolation(GridIndex<lowDimIdx> eIdx)
{
- assert(!lowDimIsBox);
+ static_assert(!isBox<lowDimIdx>(), "This interface is not available for the box method.");
lowDimElementIdx_ = eIdx;
}
- BulkPrimaryVariables interpolateBulk(const BulkSolutionVector& sol)
+ PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
{
- BulkPrimaryVariables bulkPriVars(0.0);
- if (bulkIsBox)
+ PrimaryVariables<bulkIdx> bulkPriVars(0.0);
+ if (isBox<bulkIdx>())
{
- for (std::size_t i = 0; i < bulkCornerIndices_.size(); ++i)
- for (std::size_t priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
+ for (int i = 0; i < bulkCornerIndices_.size(); ++i)
+ for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
bulkPriVars[priVarIdx] += sol[bulkCornerIndices_[i]][priVarIdx]*bulkShapeValues_[i];
}
else
@@ -106,13 +107,13 @@ public:
return bulkPriVars;
}
- LowDimPrimaryVariables interpolateLowDim(const LowDimSolutionVector& sol)
+ PrimaryVariables<lowDimIdx> interpolateLowDim(const SolutionVector<lowDimIdx>& sol) const
{
- LowDimPrimaryVariables lowDimPriVars(0.0);
- if (lowDimIsBox)
+ PrimaryVariables<lowDimIdx> lowDimPriVars(0.0);
+ if (isBox<lowDimIdx>())
{
- for (std::size_t i = 0; i < lowDimCornerIndices_.size(); ++i)
- for (std::size_t priVarIdx = 0; priVarIdx < lowDimPriVars.size(); ++priVarIdx)
+ for (int i = 0; i < lowDimCornerIndices_.size(); ++i)
+ for (int priVarIdx = 0; priVarIdx < lowDimPriVars.size(); ++priVarIdx)
lowDimPriVars[priVarIdx] += sol[lowDimCornerIndices_[i]][priVarIdx]*lowDimShapeValues_[i];
}
else
@@ -122,17 +123,18 @@ public:
return lowDimPriVars;
}
- std::size_t lowDimElementIdx() const
+ GridIndex<lowDimIdx> lowDimElementIdx() const
{ return lowDimElementIdx_; }
- std::size_t bulkElementIdx() const
+ GridIndex<bulkIdx> bulkElementIdx() const
{ return bulkElementIdx_; }
private:
ShapeValues bulkShapeValues_, lowDimShapeValues_;
- std::vector<std::size_t> bulkCornerIndices_, lowDimCornerIndices_;
- std::size_t lowDimElementIdx_;
- std::size_t bulkElementIdx_;
+ std::vector<GridIndex<bulkIdx>> bulkCornerIndices_;
+ std::vector<GridIndex<lowDimIdx>> lowDimCornerIndices_;
+ GridIndex<bulkIdx> bulkElementIdx_;
+ GridIndex<lowDimIdx> lowDimElementIdx_;
};
/*!
@@ -151,48 +153,51 @@ class PointSourceDataCircleAverage : public PointSourceData<MDTraits>
using Scalar = typename MDTraits::Scalar;
using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
- // obtain the type tags of the sub problems
- using BulkTypeTag = typename MDTraits::template SubDomain<0>::TypeTag;
- using LowDimTypeTag = typename MDTraits::template SubDomain<1>::TypeTag;
-
- using BulkPrimaryVariables = GetPropType<BulkTypeTag, Properties::PrimaryVariables>;
- using LowDimPrimaryVariables = GetPropType<LowDimTypeTag, Properties::PrimaryVariables>;
+ template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
+ template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
+ template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
+ template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
+ template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
+ template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
- using BulkSolutionVector = GetPropType<BulkTypeTag, Properties::SolutionVector>;
- using LowDimSolutionVector = GetPropType<LowDimTypeTag, Properties::SolutionVector>;
+ static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
+ static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
- enum {
- bulkIsBox = GetPropType<BulkTypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box,
- lowDimIsBox = GetPropType<LowDimTypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box
- };
+ template<std::size_t id>
+ static constexpr bool isBox()
+ { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
public:
PointSourceDataCircleAverage() : enableBulkCircleInterpolation_(false) {}
- BulkPrimaryVariables interpolateBulk(const BulkSolutionVector& sol)
+ PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
{
// bulk interpolation on the circle is only enabled for source in the
// lower dimensional domain if we use a circle distributed bulk sources
// (see coupling manager circle sources)
- BulkPrimaryVariables bulkPriVars(sol[0]);
+ PrimaryVariables<bulkIdx> bulkPriVars(sol[0]);
bulkPriVars = 0.0;
if (enableBulkCircleInterpolation_)
{
// compute the average of the bulk privars over the circle around the integration point
// this computes $\bar{p} = \frac{1}{2\pi R} int_0^2*\pi p R \text{d}\theta.
- if (bulkIsBox)
+ if (isBox<bulkIdx>())
{
assert(circleCornerIndices_.size() == circleShapeValues_.size());
Scalar weightSum = 0.0;
for (std::size_t j = 0; j < circleStencil_.size(); ++j)
{
- BulkPrimaryVariables priVars(0.0);
- const auto& cornerIndices = circleCornerIndices_[circleStencil_[j]];
- const auto& shapeValues = circleShapeValues_[circleStencil_[j]];
- for (std::size_t i = 0; i < cornerIndices.size(); ++i)
- for (std::size_t priVarIdx = 0; priVarIdx < priVars.size(); ++priVarIdx)
- priVars[priVarIdx] += sol[cornerIndices[i]][priVarIdx]*shapeValues[i];
+ PrimaryVariables<bulkIdx> priVars(0.0);
+ const auto& cornerIndices = *(circleCornerIndices_[j]);
+ const auto& shapeValues = circleShapeValues_[j];
+ for (int i = 0; i < cornerIndices.size(); ++i)
+ {
+ const auto& localSol = sol[cornerIndices[i]];
+ const auto& shapeValue = shapeValues[i];
+ for (int priVarIdx = 0; priVarIdx < PrimaryVariables<bulkIdx>::size(); ++priVarIdx)
+ priVars[priVarIdx] += localSol[priVarIdx]*shapeValue;
+ }
// multiply with weight and add
priVars *= circleIpWeight_[j];
weightSum += circleIpWeight_[j];
@@ -203,9 +208,9 @@ public:
else
{
Scalar weightSum = 0.0;
- for (std::size_t j = 0; j < circleStencil_.size(); ++j)
+ for (int j = 0; j < circleStencil_.size(); ++j)
{
- for (std::size_t priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
+ for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
bulkPriVars[priVarIdx] += sol[circleStencil_[j]][priVarIdx]*circleIpWeight_[j];
weightSum += circleIpWeight_[j];
@@ -220,10 +225,10 @@ public:
}
}
- void addCircleInterpolation(const std::unordered_map<std::size_t, std::vector<std::size_t> >& circleCornerIndices,
- const std::unordered_map<std::size_t, ShapeValues>& circleShapeValues,
+ void addCircleInterpolation(const std::vector<const std::vector<GridIndex<bulkIdx>>*>& circleCornerIndices,
+ const std::vector<ShapeValues>& circleShapeValues,
const std::vector<Scalar>& circleIpWeight,
- const std::vector<std::size_t>& circleStencil)
+ const std::vector<GridIndex<bulkIdx>>& circleStencil)
{
circleCornerIndices_ = circleCornerIndices;
circleShapeValues_ = circleShapeValues;
@@ -234,23 +239,21 @@ public:
}
void addCircleInterpolation(const std::vector<Scalar>& circleIpWeight,
- const std::vector<std::size_t>& circleStencil)
+ const std::vector<GridIndex<bulkIdx>>& circleStencil)
{
circleIpWeight_ = circleIpWeight;
circleStencil_ = circleStencil;
enableBulkCircleInterpolation_ = true;
}
- const std::vector<std::size_t>& circleStencil()
- {
- return circleStencil_;
- }
+ const std::vector<GridIndex<bulkIdx>>& circleStencil() const
+ { return circleStencil_; }
private:
- std::unordered_map<std::size_t, std::vector<std::size_t> > circleCornerIndices_;
- std::unordered_map<std::size_t, ShapeValues> circleShapeValues_;
+ std::vector<const std::vector<GridIndex<bulkIdx>>*> circleCornerIndices_;
+ std::vector<ShapeValues> circleShapeValues_;
std::vector<Scalar> circleIpWeight_;
- std::vector<std::size_t> circleStencil_;
+ std::vector<GridIndex<bulkIdx>> circleStencil_;
bool enableBulkCircleInterpolation_;
};