3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
boundary/stokesdarcy/couplingmapper.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_STOKES_DARCY_COUPLINGMAPPER_HH
26#define DUMUX_STOKES_DARCY_COUPLINGMAPPER_HH
27
28#include <type_traits>
29#include <unordered_map>
30#include <vector>
31
32#include <dune/common/exceptions.hh>
34
35namespace Dumux {
36
42{
43 struct ElementMapInfo
44 {
45 std::size_t eIdx;
46 std::size_t scvfIdx;
47 std::size_t flipScvfIdx;
48 };
49
50public:
51
55 template<class CouplingManager, class Stencils>
57 Stencils& darcyToStokesCellCenterStencils,
58 Stencils& darcyToStokesFaceStencils,
59 Stencils& stokesCellCenterToDarcyStencils,
60 Stencils& stokesFaceToDarcyStencils)
61 {
62 const auto& stokesFvGridGeometry = couplingManager.problem(CouplingManager::stokesIdx).gridGeometry();
63 const auto& darcyFvGridGeometry = couplingManager.problem(CouplingManager::darcyIdx).gridGeometry();
64
65 static_assert(std::decay_t<decltype(stokesFvGridGeometry)>::discMethod == DiscretizationMethod::staggered,
66 "The free flow domain must use the staggered discretization");
67
68 static_assert(std::decay_t<decltype(darcyFvGridGeometry)>::discMethod == DiscretizationMethod::cctpfa,
69 "The Darcy domain must use the CCTpfa discretization");
70
71 isCoupledDarcyScvf_.resize(darcyFvGridGeometry.numScvf(), false);
72
73 auto darcyFvGeometry = localView(darcyFvGridGeometry);
74 auto stokesFvGeometry = localView(stokesFvGridGeometry);
75 const auto& stokesGridView = stokesFvGridGeometry.gridView();
76
77 for (const auto& stokesElement : elements(stokesGridView))
78 {
79 stokesFvGeometry.bindElement(stokesElement);
80
81 for (const auto& scvf : scvfs(stokesFvGeometry))
82 {
83 // skip the DOF if it is not on the boundary
84 if (!scvf.boundary())
85 continue;
86
87 // get element intersecting with the scvf center
88 // for robustness, add epsilon in unit outer normal direction
89 const auto eps = (scvf.center() - stokesElement.geometry().center()).two_norm()*1e-8;
90 auto globalPos = scvf.center(); globalPos.axpy(eps, scvf.unitOuterNormal());
91 const auto darcyElementIdx = intersectingEntities(globalPos, darcyFvGridGeometry.boundingBoxTree());
92
93 // skip if no intersection was found
94 if (darcyElementIdx.empty())
95 continue;
96
97 // sanity check
98 if (darcyElementIdx.size() > 1)
99 DUNE_THROW(Dune::InvalidStateException, "Stokes face dof should only intersect with one Darcy element");
100
101 const auto stokesElementIdx = stokesFvGridGeometry.elementMapper().index(stokesElement);
102
103 const auto darcyDofIdx = darcyElementIdx[0];
104
105 stokesFaceToDarcyStencils[scvf.dofIndex()].push_back(darcyDofIdx);
106 stokesCellCenterToDarcyStencils[stokesElementIdx].push_back(darcyDofIdx);
107
108 darcyToStokesFaceStencils[darcyElementIdx[0]].push_back(scvf.dofIndex());
109 darcyToStokesCellCenterStencils[darcyElementIdx[0]].push_back(stokesElementIdx);
110
111 const auto& darcyElement = darcyFvGridGeometry.element(darcyElementIdx[0]);
112 darcyFvGeometry.bindElement(darcyElement);
113
114 // find the corresponding Darcy sub control volume face
115 for (const auto& darcyScvf : scvfs(darcyFvGeometry))
116 {
117 const auto distance = (darcyScvf.center() - scvf.center()).two_norm();
118
119 if (distance < eps)
120 {
121 isCoupledDarcyScvf_[darcyScvf.index()] = true;
122 darcyElementToStokesElementMap_[darcyElementIdx[0]].push_back({stokesElementIdx, scvf.index(), darcyScvf.index()});
123 stokesElementToDarcyElementMap_[stokesElementIdx].push_back({darcyElementIdx[0], darcyScvf.index(), scvf.index()});
124 }
125 }
126 }
127 }
128 }
129
133 bool isCoupledDarcyScvf(std::size_t darcyScvfIdx) const
134 {
135 return isCoupledDarcyScvf_[darcyScvfIdx];
136 }
137
142 {
143 return darcyElementToStokesElementMap_;
144 }
145
150 {
151 return stokesElementToDarcyElementMap_;
152 }
153
154private:
155 std::unordered_map<std::size_t, std::vector<ElementMapInfo>> darcyElementToStokesElementMap_;
156 std::unordered_map<std::size_t, std::vector<ElementMapInfo>> stokesElementToDarcyElementMap_;
157
158 std::vector<bool> isCoupledDarcyScvf_;
159};
160
161} // end namespace Dumux
162
163#endif
The available discretization methods in Dumux.
std::vector< std::size_t > intersectingEntities(const Dune::FieldVector< ctype, dimworld > &point, const BoundingBoxTree< EntitySet > &tree, bool isCartesianGrid=false)
Compute all intersections between entities and a point.
Definition: geometry/intersectingentities.hh:100
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Definition: adapt.hh:29
Coupling mapper for Stokes and Darcy domains with equal dimension.
Definition: boundary/stokesdarcy/couplingmapper.hh:42
bool isCoupledDarcyScvf(std::size_t darcyScvfIdx) const
Returns whether a Darcy scvf is coupled to the other domain.
Definition: boundary/stokesdarcy/couplingmapper.hh:133
void computeCouplingMapsAndStencils(const CouplingManager &couplingManager, Stencils &darcyToStokesCellCenterStencils, Stencils &darcyToStokesFaceStencils, Stencils &stokesCellCenterToDarcyStencils, Stencils &stokesFaceToDarcyStencils)
Main update routine.
Definition: boundary/stokesdarcy/couplingmapper.hh:56
const auto & darcyElementToStokesElementMap() const
A map that returns all Stokes elements coupled to a Darcy element.
Definition: boundary/stokesdarcy/couplingmapper.hh:141
const auto & stokesElementToDarcyElementMap() const
A map that returns all Darcy elements coupled to a Stokes element.
Definition: boundary/stokesdarcy/couplingmapper.hh:149
Definition: multidomain/couplingmanager.hh:46
const Problem< i > & problem(Dune::index_constant< i > domainIdx) const
Return a reference to the sub problem.
Definition: multidomain/couplingmanager.hh:264