26#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH
27#define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH
33#include <unordered_map>
35#include <dune/common/timer.hh>
36#include <dune/geometry/quadraturerules.hh>
50template<
class MDTraits>
54 template<std::
size_t i>
using SubDomainTypeTag =
typename MDTraits::template SubDomain<i>::TypeTag;
59 template<std::
size_t i>
63 template<std::
size_t i>
75template<
class MDTraits,
class Implementation,
class PSTraits = DefaultPo
intSourceTraits<MDTraits>>
80 using Scalar =
typename MDTraits::Scalar;
81 static constexpr auto bulkIdx =
typename MDTraits::template SubDomain<0>::Index();
82 static constexpr auto lowDimIdx =
typename MDTraits::template SubDomain<1>::Index();
83 using SolutionVector =
typename MDTraits::SolutionVector;
84 using PointSourceData =
typename PSTraits::PointSourceData;
87 template<std::
size_t id>
using PointSource =
typename PSTraits::template PointSource<id>;
88 template<std::
size_t id>
using SubDomainTypeTag =
typename MDTraits::template SubDomain<id>::TypeTag;
92 template<std::
size_t id>
using GridView =
typename GridGeometry<id>::GridView;
93 template<std::
size_t id>
using ElementMapper =
typename GridGeometry<id>::ElementMapper;
94 template<std::
size_t id>
using Element =
typename GridView<id>::template Codim<0>::Entity;
96 template<std::
size_t id>
using CouplingStencil = std::vector<GridIndex<id>>;
98 static constexpr int bulkDim = GridView<bulkIdx>::dimension;
99 static constexpr int lowDimDim = GridView<lowDimIdx>::dimension;
100 static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld;
102 template<std::
size_t id>
103 static constexpr bool isBox()
106 using GlobalPosition =
typename Element<bulkIdx>::Geometry::GlobalCoordinate;
115 template<std::
size_t id>
using CouplingStencils = std::unordered_map<GridIndex<id>, CouplingStencil<id>>;
121 std::shared_ptr<
const GridGeometry<lowDimIdx>> lowDimGridGeometry)
123 glue_ = std::make_shared<GlueType>();
130 std::shared_ptr<
const GridGeometry<lowDimIdx>> lowDimGridGeometry)
140 void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem,
141 std::shared_ptr<Problem<lowDimIdx>> lowDimProblem,
142 const SolutionVector&
curSol)
145 this->
setSubProblems(std::make_tuple(bulkProblem, lowDimProblem));
147 integrationOrder_ = getParam<int>(
"MixedDimension.IntegrationOrder", 1);
148 asImp_().computePointSourceData(integrationOrder_);
172 template<std::
size_t i, std::
size_t j>
174 const Element<i>& element,
175 Dune::index_constant<j> domainJ)
const
177 static_assert(i != j,
"A domain cannot be coupled to itself!");
179 const auto eIdx = this->
problem(domainI).gridGeometry().elementMapper().index(element);
198 template<std::
size_t i, std::
size_t j,
class LocalAssemblerI>
200 const LocalAssemblerI& localAssemblerI,
201 Dune::index_constant<j> domainJ,
202 std::size_t dofIdxGlobalJ)
204 static_assert(i != j,
"A domain cannot be coupled to itself!");
206 typename LocalAssemblerI::LocalResidual::ElementResidualVector residual;
208 const auto& element = localAssemblerI.element();
209 const auto& fvGeometry = localAssemblerI.fvGeometry();
210 const auto& curElemVolVars = localAssemblerI.curElemVolVars();
212 residual.resize(fvGeometry.numScv());
213 for (
const auto& scv : scvs(fvGeometry))
215 auto couplingSource = this->
problem(domainI).scvPointSources(element, fvGeometry, curElemVolVars, scv);
216 couplingSource += this->
problem(domainI).source(element, fvGeometry, curElemVolVars, scv);
217 couplingSource *= -GridGeometry<i>::Extrusion::volume(scv)*curElemVolVars[scv].extrusionFactor();
218 residual[scv.indexInElement()] = couplingSource;
239 std::cout <<
"Initializing the point sources..." << std::endl;
249 const auto& bulkGridGeometry = this->
problem(bulkIdx).gridGeometry();
250 const auto& lowDimGridGeometry = this->
problem(lowDimIdx).gridGeometry();
255 pointSourceData_.reserve(this->
glue().size());
256 averageDistanceToBulkCell_.reserve(this->
glue().size());
257 for (
const auto& is : intersections(this->
glue()))
260 const auto& inside = is.targetEntity(0);
262 const auto intersectionGeometry = is.geometry();
265 const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
266 const std::size_t lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
269 for (
auto&& qp : quad)
272 for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
274 const auto& outside = is.domainEntity(outsideIdx);
275 const std::size_t bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
278 const auto globalPos = intersectionGeometry.global(qp.position());
280 const auto qpweight = qp.weight();
281 const auto ie = intersectionGeometry.integrationElement(qp.position());
282 pointSources(bulkIdx).emplace_back(globalPos,
id, qpweight, ie, std::vector<std::size_t>({bulkElementIdx}));
283 pointSources(bulkIdx).back().setEmbeddings(is.numDomainNeighbors());
284 pointSources(lowDimIdx).emplace_back(globalPos,
id, qpweight, ie, std::vector<std::size_t>({lowDimElementIdx}));
285 pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors());
289 PointSourceData psData;
291 if constexpr (isBox<lowDimIdx>())
293 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
294 const auto lowDimGeometry = this->
problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
295 ShapeValues shapeValues;
296 this->
getShapeValues(lowDimIdx, this->
problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
297 psData.addLowDimInterpolation(shapeValues, this->
vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
301 psData.addLowDimInterpolation(lowDimElementIdx);
305 if constexpr (isBox<bulkIdx>())
307 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
308 const auto bulkGeometry = this->
problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
309 ShapeValues shapeValues;
310 this->
getShapeValues(bulkIdx, this->
problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues);
311 psData.addBulkInterpolation(shapeValues, this->
vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
315 psData.addBulkInterpolation(bulkElementIdx);
326 if (isBox<bulkIdx>())
328 const auto& vertices = this->
vertexIndices(bulkIdx, bulkElementIdx);
330 vertices.begin(), vertices.end());
334 this->
couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx);
339 if (isBox<lowDimIdx>())
341 const auto& vertices = this->
vertexIndices(lowDimIdx, lowDimElementIdx);
343 vertices.begin(), vertices.end());
348 this->
couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx);
355 using namespace Dune::Hybrid;
356 forEach(integralRange(std::integral_constant<std::size_t, MDTraits::numSubDomains>{}), [&](
const auto domainIdx)
360 std::sort(stencil.second.begin(), stencil.second.end());
361 stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end());
365 std::cout <<
"took " << watch.elapsed() <<
" seconds." << std::endl;
375 {
return pointSourceData_[id]; }
378 template<std::
size_t id>
379 const GridView<id>&
gridView(Dune::index_constant<id> domainIdx)
const
380 {
return this->
problem(domainIdx).gridGeometry().gridView(); }
384 {
return pointSourceData_[id].interpolateBulk(this->
curSol()[bulkIdx]); }
388 {
return pointSourceData_[id].interpolateLowDim(this->
curSol()[lowDimIdx]); }
392 {
return averageDistanceToBulkCell_[id]; }
396 {
return std::get<bulkIdx>(pointSources_); }
400 {
return std::get<lowDimIdx>(pointSources_); }
403 template<std::
size_t i>
404 const std::vector<PointSource<i>>&
pointSources(Dune::index_constant<i> dom)
const
405 {
return std::get<i>(pointSources_); }
408 template<std::
size_t i>
410 {
return std::get<i>(couplingStencils_); }
414 {
return pointSourceData_; }
417 template<std::
size_t i>
418 const CouplingStencil<i>&
emptyStencil(Dune::index_constant<i> dom)
const
419 {
return std::get<i>(emptyStencil_); }
424 template<std::
size_t id>
428 if constexpr (isBox<domainIdx>())
431 for (
const auto& element : elements(
gridView(domainIdx)))
433 constexpr int dim = GridView<domainIdx>::dimension;
434 const auto eIdx = this->
problem(domainIdx).gridGeometry().elementMapper().index(element);
435 this->
vertexIndices(domainIdx, eIdx).resize(element.subEntities(dim));
436 for (
int i = 0; i < element.subEntities(dim); ++i)
437 this->
vertexIndices(domainIdx, eIdx)[i] = this->
problem(domainIdx).gridGeometry().vertexMapper().subIndex(element, i, dim);
443 template<std::
size_t i,
class FVGG,
class Geometry,
class ShapeValues>
444 void getShapeValues(Dune::index_constant<i> domainI,
const FVGG& gridGeometry,
const Geometry& geo,
const GlobalPosition& globalPos, ShapeValues& shapeValues)
448 const auto ipLocal = geo.local(globalPos);
449 const auto& localBasis = this->
problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis();
450 localBasis.evaluateFunction(ipLocal, shapeValues);
453 DUNE_THROW(Dune::InvalidStateException,
"Shape values requested for other discretization than box!");
465 pointSourceData_.clear();
466 averageDistanceToBulkCell_.clear();
474 const auto& bulkGridGeometry = this->
problem(bulkIdx).gridGeometry();
475 const auto& lowDimGridGeometry = this->
problem(lowDimIdx).gridGeometry();
478 glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree());
483 {
return pointSourceData_; }
487 {
return averageDistanceToBulkCell_; }
490 template<std::
size_t i>
492 {
return std::get<i>(pointSources_); }
495 template<std::
size_t i>
497 {
return std::get<i>(couplingStencils_); }
500 template<std::
size_t i>
501 std::vector<GridIndex<i>>&
vertexIndices(Dune::index_constant<i> dom, GridIndex<i> eIdx)
502 {
return std::get<i>(vertexIndices_)[eIdx]; }
505 template<std::
size_t i>
506 std::vector<std::vector<GridIndex<i>>>&
vertexIndices(Dune::index_constant<i> dom)
507 {
return std::get<i>(vertexIndices_); }
514 {
return *
static_cast<Implementation *
>(
this); }
518 {
return *
static_cast<const Implementation *
>(
this); }
526 std::tuple<std::vector<PointSource<bulkIdx>>, std::vector<PointSource<lowDimIdx>>> pointSources_;
527 std::vector<PointSourceData> pointSourceData_;
528 std::vector<Scalar> averageDistanceToBulkCell_;
531 std::tuple<std::vector<std::vector<GridIndex<bulkIdx>>>,
532 std::vector<std::vector<GridIndex<lowDimIdx>>>> vertexIndices_;
534 std::tuple<CouplingStencil<bulkIdx>, CouplingStencil<lowDimIdx>> emptyStencil_;
537 std::shared_ptr<GlueType> glue_;
540 int integrationOrder_ = 1;
The available discretization methods in Dumux.
An integration point source class, i.e. sources located at a single point in space associated with a ...
Data associated with a point source.
A class glueing two grids of potentially different dimension geometrically. Intersections are compute...
Geometry::ctype averageDistancePointGeometry(const typename Geometry::GlobalCoordinate &p, const Geometry &geometry, std::size_t integrationOrder=2)
Compute the average distance from a point to a geometry by integration.
Definition: geometry/distance.hh:36
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Struture to define the index types used for grid and local indices.
Definition: indextraits.hh:38
A vector of primary variables.
Definition: common/properties.hh:49
A vector of size number equations that can be used for Neumann fluxes, sources, residuals,...
Definition: common/properties.hh:51
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:57
Definition: common/properties.hh:101
A class representing the intersection entites two geometric entity sets.
Definition: geometry/intersectionentityset.hh:55
Definition: multidomain/couplingmanager.hh:46
void setSubProblems(const std::tuple< std::shared_ptr< SubProblems >... > &problems)
set the pointers to the sub problems
Definition: multidomain/couplingmanager.hh:247
const Problem< i > & problem(Dune::index_constant< i > domainIdx) const
Return a reference to the sub problem.
Definition: multidomain/couplingmanager.hh:264
SolutionVector & curSol()
the solution vector of the coupled problem
Definition: multidomain/couplingmanager.hh:278
void updateSolution(const SolutionVector &curSol)
Updates the entire solution vector, e.g. before assembly or after grid adaption.
Definition: multidomain/couplingmanager.hh:186
the default point source traits
Definition: couplingmanagerbase.hh:52
Manages the coupling between bulk elements and lower dimensional elements Point sources on each integ...
Definition: couplingmanagerbase.hh:78
std::vector< GridIndex< i > > & vertexIndices(Dune::index_constant< i > dom, GridIndex< i > eIdx)
Return a reference to the vertex indices.
Definition: couplingmanagerbase.hh:501
Scalar averageDistance(std::size_t id) const
return the average distance to the coupled bulk cell center
Definition: couplingmanagerbase.hh:391
Implementation & asImp_()
Returns the implementation of the problem (i.e. static polymorphism)
Definition: couplingmanagerbase.hh:513
const std::vector< PointSource< i > > & pointSources(Dune::index_constant< i > dom) const
Return the point source if domain i.
Definition: couplingmanagerbase.hh:404
PSTraits PointSourceTraits
export the point source traits
Definition: couplingmanagerbase.hh:113
EmbeddedCouplingManagerBase(std::shared_ptr< const GridGeometry< bulkIdx > > bulkGridGeometry, std::shared_ptr< const GridGeometry< lowDimIdx > > lowDimGridGeometry)
Constructor.
Definition: couplingmanagerbase.hh:129
const CouplingStencils< i > & couplingStencils(Dune::index_constant< i > dom) const
Return reference to bulk coupling stencil member of domain i.
Definition: couplingmanagerbase.hh:409
void updateAfterGridAdaption(std::shared_ptr< const GridGeometry< bulkIdx > > bulkGridGeometry, std::shared_ptr< const GridGeometry< lowDimIdx > > lowDimGridGeometry)
call this after grid adaption
Definition: couplingmanagerbase.hh:120
const GridView< id > & gridView(Dune::index_constant< id > domainIdx) const
Return a reference to the bulk problem.
Definition: couplingmanagerbase.hh:379
const Implementation & asImp_() const
Returns the implementation of the problem (i.e. static polymorphism)
Definition: couplingmanagerbase.hh:517
std::vector< Scalar > & averageDistanceToBulkCell()
Return reference to average distances to bulk cell.
Definition: couplingmanagerbase.hh:486
const CouplingStencil< i > & emptyStencil(Dune::index_constant< i > dom) const
Return a reference to an empty stencil.
Definition: couplingmanagerbase.hh:418
const std::vector< PointSource< bulkIdx > > & bulkPointSources() const
Return reference to bulk point sources.
Definition: couplingmanagerbase.hh:395
const PointSourceData & pointSourceData(std::size_t id) const
Methods to be accessed by the subproblems.
Definition: couplingmanagerbase.hh:374
std::vector< std::vector< GridIndex< i > > > & vertexIndices(Dune::index_constant< i > dom)
Return a reference to the vertex indices container.
Definition: couplingmanagerbase.hh:506
void glueGrids()
compute the intersections between the two grids
Definition: couplingmanagerbase.hh:472
MDTraits MultiDomainTraits
export traits
Definition: couplingmanagerbase.hh:111
void computePointSourceData(std::size_t order=1, bool verbose=false)
Definition: couplingmanagerbase.hh:234
std::vector< PointSourceData > & pointSourceData()
Return reference to point source data vector member.
Definition: couplingmanagerbase.hh:482
const CouplingStencil< j > & couplingStencil(Dune::index_constant< i > domainI, const Element< i > &element, Dune::index_constant< j > domainJ) const
Methods to be accessed by the assembly.
Definition: couplingmanagerbase.hh:173
PrimaryVariables< lowDimIdx > lowDimPriVars(std::size_t id) const
Return data for a low dim point source with the identifier id.
Definition: couplingmanagerbase.hh:387
void getShapeValues(Dune::index_constant< i > domainI, const FVGG &gridGeometry, const Geometry &geo, const GlobalPosition &globalPos, ShapeValues &shapeValues)
compute the shape function for a given point and geometry
Definition: couplingmanagerbase.hh:444
const std::vector< PointSource< lowDimIdx > > & lowDimPointSources() const
Return reference to low dim point sources.
Definition: couplingmanagerbase.hh:399
void clear()
Clear all internal data members.
Definition: couplingmanagerbase.hh:457
const GlueType & glue() const
Definition: couplingmanagerbase.hh:509
decltype(auto) evalCouplingResidual(Dune::index_constant< i > domainI, const LocalAssemblerI &localAssemblerI, Dune::index_constant< j > domainJ, std::size_t dofIdxGlobalJ)
evaluates the element residual of a coupled element of domain i which depends on the variables at the...
Definition: couplingmanagerbase.hh:199
void init(std::shared_ptr< Problem< bulkIdx > > bulkProblem, std::shared_ptr< Problem< lowDimIdx > > lowDimProblem, const SolutionVector &curSol)
Methods to be accessed by main.
Definition: couplingmanagerbase.hh:140
CouplingStencils< i > & couplingStencils(Dune::index_constant< i > dom)
Return reference to bulk coupling stencil member of domain i.
Definition: couplingmanagerbase.hh:496
std::vector< PointSource< i > > & pointSources(Dune::index_constant< i > dom)
Return the point source if domain i.
Definition: couplingmanagerbase.hh:491
std::unordered_map< GridIndex< id >, CouplingStencil< id > > CouplingStencils
export stencil types
Definition: couplingmanagerbase.hh:115
PrimaryVariables< bulkIdx > bulkPriVars(std::size_t id) const
Return data for a bulk point source with the identifier id.
Definition: couplingmanagerbase.hh:383
void precomputeVertexIndices(Dune::index_constant< id > domainIdx)
computes the vertex indices per element for the box method
Definition: couplingmanagerbase.hh:425
std::size_t idCounter_
id generator for point sources
Definition: couplingmanagerbase.hh:521
const std::vector< PointSourceData > & pointSourceData() const
Return reference to point source data vector member.
Definition: couplingmanagerbase.hh:413
An integration point source class with an identifier to attach data and a quadrature weight and integ...
Definition: integrationpointsource.hh:45
A helper class calculating a DOF-index to point source map.
Definition: integrationpointsource.hh:135
A point source data class used for integration in multidimension models.
Definition: pointsourcedata.hh:43
Helper functions for distance queries.
Algorithms that finds which geometric entites intersect.
Declares all properties used in Dumux.
The interface of the coupling manager for multi domain problems.