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>
51template<
class MDTraits>
55 template<std::
size_t i>
using SubDomainTypeTag =
typename MDTraits::template SubDomain<i>::TypeTag;
60 template<std::
size_t i>
64 template<std::
size_t i>
76template<
class MDTraits,
class Implementation,
class PSTraits = DefaultPo
intSourceTraits<MDTraits>>
81 using Scalar =
typename MDTraits::Scalar;
82 static constexpr auto bulkIdx =
typename MDTraits::template SubDomain<0>::Index();
83 static constexpr auto lowDimIdx =
typename MDTraits::template SubDomain<1>::Index();
84 using SolutionVector =
typename MDTraits::SolutionVector;
85 using PointSourceData =
typename PSTraits::PointSourceData;
88 template<std::
size_t id>
using PointSource =
typename PSTraits::template PointSource<id>;
89 template<std::
size_t id>
using SubDomainTypeTag =
typename MDTraits::template SubDomain<id>::TypeTag;
93 template<std::
size_t id>
using GridView =
typename GridGeometry<id>::GridView;
94 template<std::
size_t id>
using ElementMapper =
typename GridGeometry<id>::ElementMapper;
95 template<std::
size_t id>
using Element =
typename GridView<id>::template Codim<0>::Entity;
97 template<std::
size_t id>
using CouplingStencil = std::vector<GridIndex<id>>;
99 static constexpr int bulkDim = GridView<bulkIdx>::dimension;
100 static constexpr int lowDimDim = GridView<lowDimIdx>::dimension;
101 static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld;
103 template<std::
size_t id>
104 static constexpr bool isBox()
107 using GlobalPosition =
typename Element<bulkIdx>::Geometry::GlobalCoordinate;
116 template<std::
size_t id>
using CouplingStencils = std::unordered_map<GridIndex<id>, CouplingStencil<id>>;
122 std::shared_ptr<
const GridGeometry<lowDimIdx>> lowDimGridGeometry)
124 glue_ = std::make_shared<GlueType>();
131 std::shared_ptr<
const GridGeometry<lowDimIdx>> lowDimGridGeometry)
141 void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem,
142 std::shared_ptr<Problem<lowDimIdx>> lowDimProblem,
143 const SolutionVector&
curSol)
146 this->
setSubProblems(std::make_tuple(bulkProblem, lowDimProblem));
148 integrationOrder_ = getParam<int>(
"MixedDimension.IntegrationOrder", 1);
149 asImp_().computePointSourceData(integrationOrder_);
173 template<std::
size_t i, std::
size_t j>
175 const Element<i>& element,
176 Dune::index_constant<j> domainJ)
const
178 static_assert(i != j,
"A domain cannot be coupled to itself!");
180 const auto eIdx = this->
problem(domainI).gridGeometry().elementMapper().index(element);
199 template<std::
size_t i, std::
size_t j,
class LocalAssemblerI>
201 const LocalAssemblerI& localAssemblerI,
202 Dune::index_constant<j> domainJ,
203 std::size_t dofIdxGlobalJ)
205 static_assert(i != j,
"A domain cannot be coupled to itself!");
207 typename LocalAssemblerI::LocalResidual::ElementResidualVector residual;
209 const auto& element = localAssemblerI.element();
210 const auto& fvGeometry = localAssemblerI.fvGeometry();
211 const auto& curElemVolVars = localAssemblerI.curElemVolVars();
213 residual.resize(fvGeometry.numScv());
214 for (
const auto& scv : scvs(fvGeometry))
216 auto couplingSource = this->
problem(domainI).scvPointSources(element, fvGeometry, curElemVolVars, scv);
217 couplingSource += this->
problem(domainI).source(element, fvGeometry, curElemVolVars, scv);
219 residual[scv.indexInElement()] = couplingSource;
240 std::cout <<
"Initializing the point sources..." << std::endl;
250 const auto& bulkGridGeometry = this->
problem(bulkIdx).gridGeometry();
251 const auto& lowDimGridGeometry = this->
problem(lowDimIdx).gridGeometry();
256 pointSourceData_.reserve(this->
glue().size());
257 averageDistanceToBulkCell_.reserve(this->
glue().size());
258 for (
const auto& is : intersections(this->
glue()))
261 const auto& inside = is.targetEntity(0);
263 const auto intersectionGeometry = is.geometry();
266 const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
267 const std::size_t lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
270 for (
auto&& qp : quad)
273 for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
275 const auto& outside = is.domainEntity(outsideIdx);
276 const std::size_t bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
279 const auto globalPos = intersectionGeometry.global(qp.position());
281 const auto qpweight = qp.weight();
282 const auto ie = intersectionGeometry.integrationElement(qp.position());
283 pointSources(bulkIdx).emplace_back(globalPos,
id, qpweight, ie, bulkElementIdx);
284 pointSources(bulkIdx).back().setEmbeddings(is.numDomainNeighbors());
285 pointSources(lowDimIdx).emplace_back(globalPos,
id, qpweight, ie, lowDimElementIdx);
286 pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors());
290 PointSourceData psData;
292 if constexpr (isBox<lowDimIdx>())
294 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
295 const auto lowDimGeometry = this->
problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
296 ShapeValues shapeValues;
297 this->
getShapeValues(lowDimIdx, this->
problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
298 psData.addLowDimInterpolation(shapeValues, this->
vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
302 psData.addLowDimInterpolation(lowDimElementIdx);
306 if constexpr (isBox<bulkIdx>())
308 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
309 const auto bulkGeometry = this->
problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
310 ShapeValues shapeValues;
311 this->
getShapeValues(bulkIdx, this->
problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues);
312 psData.addBulkInterpolation(shapeValues, this->
vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
316 psData.addBulkInterpolation(bulkElementIdx);
327 if (isBox<bulkIdx>())
329 const auto& vertices = this->
vertexIndices(bulkIdx, bulkElementIdx);
331 vertices.begin(), vertices.end());
335 this->
couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx);
340 if (isBox<lowDimIdx>())
342 const auto& vertices = this->
vertexIndices(lowDimIdx, lowDimElementIdx);
344 vertices.begin(), vertices.end());
349 this->
couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx);
356 using namespace Dune::Hybrid;
357 forEach(integralRange(std::integral_constant<std::size_t, MDTraits::numSubDomains>{}), [&](
const auto domainIdx)
361 std::sort(stencil.second.begin(), stencil.second.end());
362 stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end());
366 std::cout <<
"took " << watch.elapsed() <<
" seconds." << std::endl;
376 {
return pointSourceData_[id]; }
379 template<std::
size_t id>
380 const GridView<id>&
gridView(Dune::index_constant<id> domainIdx)
const
381 {
return this->
problem(domainIdx).gridGeometry().gridView(); }
385 {
return pointSourceData_[id].interpolateBulk(this->
curSol()[bulkIdx]); }
389 {
return pointSourceData_[id].interpolateLowDim(this->
curSol()[lowDimIdx]); }
393 {
return averageDistanceToBulkCell_[id]; }
397 {
return std::get<bulkIdx>(pointSources_); }
401 {
return std::get<lowDimIdx>(pointSources_); }
404 template<std::
size_t i>
405 const std::vector<PointSource<i>>&
pointSources(Dune::index_constant<i> dom)
const
406 {
return std::get<i>(pointSources_); }
409 template<std::
size_t i>
411 {
return std::get<i>(couplingStencils_); }
415 {
return pointSourceData_; }
418 template<std::
size_t i>
419 const CouplingStencil<i>&
emptyStencil(Dune::index_constant<i> dom)
const
420 {
return std::get<i>(emptyStencil_); }
425 template<std::
size_t id>
429 if constexpr (isBox<domainIdx>())
432 for (
const auto& element : elements(
gridView(domainIdx)))
434 constexpr int dim = GridView<domainIdx>::dimension;
435 const auto eIdx = this->
problem(domainIdx).gridGeometry().elementMapper().index(element);
436 this->
vertexIndices(domainIdx, eIdx).resize(element.subEntities(dim));
437 for (
int i = 0; i < element.subEntities(dim); ++i)
438 this->
vertexIndices(domainIdx, eIdx)[i] = this->
problem(domainIdx).gridGeometry().vertexMapper().subIndex(element, i, dim);
444 template<std::
size_t i,
class FVGG,
class Geometry,
class ShapeValues>
445 void getShapeValues(Dune::index_constant<i> domainI,
const FVGG& gridGeometry,
const Geometry& geo,
const GlobalPosition& globalPos, ShapeValues& shapeValues)
449 const auto ipLocal = geo.local(globalPos);
450 const auto& localBasis = this->
problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis();
451 localBasis.evaluateFunction(ipLocal, shapeValues);
454 DUNE_THROW(Dune::InvalidStateException,
"Shape values requested for other discretization than box!");
466 pointSourceData_.clear();
467 averageDistanceToBulkCell_.clear();
475 const auto& bulkGridGeometry = this->
problem(bulkIdx).gridGeometry();
476 const auto& lowDimGridGeometry = this->
problem(lowDimIdx).gridGeometry();
479 glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree());
484 {
return pointSourceData_; }
488 {
return averageDistanceToBulkCell_; }
491 template<std::
size_t i>
493 {
return std::get<i>(pointSources_); }
496 template<std::
size_t i>
498 {
return std::get<i>(couplingStencils_); }
501 template<std::
size_t i>
502 std::vector<GridIndex<i>>&
vertexIndices(Dune::index_constant<i> dom, GridIndex<i> eIdx)
503 {
return std::get<i>(vertexIndices_)[eIdx]; }
506 template<std::
size_t i>
507 std::vector<std::vector<GridIndex<i>>>&
vertexIndices(Dune::index_constant<i> dom)
508 {
return std::get<i>(vertexIndices_); }
515 {
return *
static_cast<Implementation *
>(
this); }
519 {
return *
static_cast<const Implementation *
>(
this); }
527 std::tuple<std::vector<PointSource<bulkIdx>>, std::vector<PointSource<lowDimIdx>>> pointSources_;
528 std::vector<PointSourceData> pointSourceData_;
529 std::vector<Scalar> averageDistanceToBulkCell_;
532 std::tuple<std::vector<std::vector<GridIndex<bulkIdx>>>,
533 std::vector<std::vector<GridIndex<lowDimIdx>>>> vertexIndices_;
535 std::tuple<CouplingStencil<bulkIdx>, CouplingStencil<lowDimIdx>> emptyStencil_;
538 std::shared_ptr<GlueType> glue_;
541 int integrationOrder_ = 1;
A helper to deduce a vector with the same size as numbers of equations.
The available discretization methods in Dumux.
Algorithms that finds which geometric entites intersect.
Helper functions for distance queries.
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: distance.hh:36
typename NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:46
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
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
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:57
Definition: common/properties.hh:98
A class representing the intersection entites two geometric entity sets.
Definition: 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:53
Manages the coupling between bulk elements and lower dimensional elements Point sources on each integ...
Definition: couplingmanagerbase.hh:79
std::vector< GridIndex< i > > & vertexIndices(Dune::index_constant< i > dom, GridIndex< i > eIdx)
Return a reference to the vertex indices.
Definition: couplingmanagerbase.hh:502
Scalar averageDistance(std::size_t id) const
return the average distance to the coupled bulk cell center
Definition: couplingmanagerbase.hh:392
Implementation & asImp_()
Returns the implementation of the problem (i.e. static polymorphism)
Definition: couplingmanagerbase.hh:514
const std::vector< PointSource< i > > & pointSources(Dune::index_constant< i > dom) const
Return the point source if domain i.
Definition: couplingmanagerbase.hh:405
PSTraits PointSourceTraits
export the point source traits
Definition: couplingmanagerbase.hh:114
EmbeddedCouplingManagerBase(std::shared_ptr< const GridGeometry< bulkIdx > > bulkGridGeometry, std::shared_ptr< const GridGeometry< lowDimIdx > > lowDimGridGeometry)
Constructor.
Definition: couplingmanagerbase.hh:130
const CouplingStencils< i > & couplingStencils(Dune::index_constant< i > dom) const
Return reference to bulk coupling stencil member of domain i.
Definition: couplingmanagerbase.hh:410
void updateAfterGridAdaption(std::shared_ptr< const GridGeometry< bulkIdx > > bulkGridGeometry, std::shared_ptr< const GridGeometry< lowDimIdx > > lowDimGridGeometry)
call this after grid adaption
Definition: couplingmanagerbase.hh:121
const GridView< id > & gridView(Dune::index_constant< id > domainIdx) const
Return a reference to the bulk problem.
Definition: couplingmanagerbase.hh:380
const Implementation & asImp_() const
Returns the implementation of the problem (i.e. static polymorphism)
Definition: couplingmanagerbase.hh:518
std::vector< Scalar > & averageDistanceToBulkCell()
Return reference to average distances to bulk cell.
Definition: couplingmanagerbase.hh:487
const CouplingStencil< i > & emptyStencil(Dune::index_constant< i > dom) const
Return a reference to an empty stencil.
Definition: couplingmanagerbase.hh:419
const std::vector< PointSource< bulkIdx > > & bulkPointSources() const
Return reference to bulk point sources.
Definition: couplingmanagerbase.hh:396
const PointSourceData & pointSourceData(std::size_t id) const
Methods to be accessed by the subproblems.
Definition: couplingmanagerbase.hh:375
std::vector< std::vector< GridIndex< i > > > & vertexIndices(Dune::index_constant< i > dom)
Return a reference to the vertex indices container.
Definition: couplingmanagerbase.hh:507
void glueGrids()
compute the intersections between the two grids
Definition: couplingmanagerbase.hh:473
MDTraits MultiDomainTraits
export traits
Definition: couplingmanagerbase.hh:112
void computePointSourceData(std::size_t order=1, bool verbose=false)
Definition: couplingmanagerbase.hh:235
std::vector< PointSourceData > & pointSourceData()
Return reference to point source data vector member.
Definition: couplingmanagerbase.hh:483
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:174
PrimaryVariables< lowDimIdx > lowDimPriVars(std::size_t id) const
Return data for a low dim point source with the identifier id.
Definition: couplingmanagerbase.hh:388
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:445
const std::vector< PointSource< lowDimIdx > > & lowDimPointSources() const
Return reference to low dim point sources.
Definition: couplingmanagerbase.hh:400
void clear()
Clear all internal data members.
Definition: couplingmanagerbase.hh:458
const GlueType & glue() const
Definition: couplingmanagerbase.hh:510
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:200
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:141
CouplingStencils< i > & couplingStencils(Dune::index_constant< i > dom)
Return reference to bulk coupling stencil member of domain i.
Definition: couplingmanagerbase.hh:497
std::vector< PointSource< i > > & pointSources(Dune::index_constant< i > dom)
Return the point source if domain i.
Definition: couplingmanagerbase.hh:492
std::unordered_map< GridIndex< id >, CouplingStencil< id > > CouplingStencils
export stencil types
Definition: couplingmanagerbase.hh:116
PrimaryVariables< bulkIdx > bulkPriVars(std::size_t id) const
Return data for a bulk point source with the identifier id.
Definition: couplingmanagerbase.hh:384
void precomputeVertexIndices(Dune::index_constant< id > domainIdx)
computes the vertex indices per element for the box method
Definition: couplingmanagerbase.hh:426
std::size_t idCounter_
id generator for point sources
Definition: couplingmanagerbase.hh:522
const std::vector< PointSourceData > & pointSourceData() const
Return reference to point source data vector member.
Definition: couplingmanagerbase.hh:414
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:112
A point source data class used for integration in multidimension models.
Definition: pointsourcedata.hh:43
Declares all properties used in Dumux.
The interface of the coupling manager for multi domain problems.