3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
boundary/darcydarcy/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_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
26#define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
27
28#include <iostream>
29#include <unordered_map>
30#include <tuple>
31#include <vector>
32
33#include <dune/common/timer.hh>
34#include <dune/common/exceptions.hh>
35#include <dune/common/indices.hh>
36
39
40namespace Dumux {
41
47template<class MDTraits>
49{
50 using Scalar = typename MDTraits::Scalar;
51
52 template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry;
53 template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace;
54 template<std::size_t i> using GridView = typename GridGeometry<i>::GridView;
55 template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
56
57 template<std::size_t i>
58 static constexpr auto domainIdx()
59 { return typename MDTraits::template SubDomain<i>::Index{}; }
60
61 template<std::size_t i>
62 static constexpr bool isCCTpfa()
63 { return GridGeometry<i>::discMethod == DiscretizationMethod::cctpfa; }
64
65 struct ScvfInfo
66 {
67 std::size_t eIdxOutside;
68 std::size_t flipScvfIdx;
69 };
70
71 using FlipScvfMapType = std::unordered_map<std::size_t, ScvfInfo>;
72 using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>;
73
74 static constexpr std::size_t numSD = MDTraits::numSubDomains;
75
76public:
80 template<class CouplingManager>
81 void update(const CouplingManager& couplingManager)
82 {
83 // TODO: Box and multiple domains
84 static_assert(numSD == 2, "More than two subdomains not implemented!");
85 static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!");
86
87 Dune::Timer watch;
88 std::cout << "Initializing the coupling map..." << std::endl;
89
90 for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
91 {
92 stencils_[domIdx].clear();
93 scvfInfo_[domIdx].clear();
94 }
95
96 const auto& problem0 = couplingManager.problem(domainIdx<0>());
97 const auto& problem1 = couplingManager.problem(domainIdx<1>());
98 const auto& gg0 = problem0.gridGeometry();
99 const auto& gg1 = problem1.gridGeometry();
100
101 isCoupledScvf_[0].resize(gg0.numScvf(), false);
102 isCoupledScvf_[1].resize(gg1.numScvf(), false);
103
104 for (const auto& element0 : elements(gg0.gridView()))
105 {
106 auto fvGeometry0 = localView(gg0);
107 fvGeometry0.bindElement(element0);
108
109 for (const auto& scvf0 : scvfs(fvGeometry0))
110 {
111 // skip all non-boundaries
112 if (!scvf0.boundary())
113 continue;
114
115 // get elements intersecting with the scvf center
116 // for robustness add epsilon in unit outer normal direction
117 const auto eps = (scvf0.ipGlobal() - element0.geometry().corner(0)).two_norm()*1e-8;
118 auto globalPos = scvf0.ipGlobal(); globalPos.axpy(eps, scvf0.unitOuterNormal());
119 const auto indices = intersectingEntities(globalPos, gg1.boundingBoxTree());
120
121 // skip if no intersection was found
122 if (indices.empty())
123 continue;
124
125 // sanity check
126 if (indices.size() > 1)
127 DUNE_THROW(Dune::InvalidStateException, "Are you sure your grids is conforming at the boundary?");
128
129 // add the pair to the multimap
130 const auto eIdx0 = gg0.elementMapper().index(element0);
131 const auto eIdx1 = indices[0];
132 stencils_[0][eIdx0].push_back(eIdx1);
133
134 // mark the scvf and find and mark the flip scvf
135 isCoupledScvf_[0][scvf0.index()] = true;
136 const auto& element1 = gg1.element(eIdx1);
137 auto fvGeometry1 = localView(gg1);
138 fvGeometry1.bindElement(element1);
139
140 using std::abs;
141 for (const auto& scvf1 : scvfs(fvGeometry1))
142 if (scvf1.boundary())
143 if (abs(scvf1.unitOuterNormal()*scvf0.unitOuterNormal() + 1) < eps)
144 {
145 isCoupledScvf_[1][scvf1.index()] = true;
146 scvfInfo_[0][scvf0.index()] = ScvfInfo{eIdx1, scvf1.index()};
147 scvfInfo_[1][scvf1.index()] = ScvfInfo{eIdx0, scvf0.index()};
148 }
149 }
150 }
151
152 // create the inverse map for efficient access
153 for (const auto& entry : stencils_[0])
154 for (const auto idx : entry.second)
155 stencils_[1][idx].push_back(entry.first);
156
157 std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
158 }
159
174 template<std::size_t i, std::size_t j>
175 const std::vector<std::size_t>& couplingStencil(Dune::index_constant<i> domainI,
176 const std::size_t eIdxI,
177 Dune::index_constant<j> domainJ) const
178 {
179 static_assert(i != j, "A domain cannot be coupled to itself!");
180 if (isCoupledElement(domainI, eIdxI))
181 return stencils_[i].at(eIdxI);
182 else
183 return emptyStencil_;
184 }
185
189 template<std::size_t i>
190 bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
191 { return static_cast<bool>(stencils_[i].count(eIdx)); }
192
198 template<std::size_t i>
199 bool isCoupled(Dune::index_constant<i> domainI,
200 const SubControlVolumeFace<i>& scvf) const
201 {
202 return isCoupledScvf_[i].at(scvf.index());
203 }
204
210 template<std::size_t i>
211 std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
212 const SubControlVolumeFace<i>& scvf) const
213 {
214 return scvfInfo_[i].at(scvf.index()).flipScvfIdx;
215 }
216
222 template<std::size_t i>
223 std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
224 const SubControlVolumeFace<i>& scvf) const
225 {
226 return scvfInfo_[i].at(scvf.index()).eIdxOutside;
227 }
228
229private:
230 std::array<MapType, numSD> stencils_;
231 std::vector<std::size_t> emptyStencil_;
232 std::array<std::vector<bool>, numSD> isCoupledScvf_;
233 std::array<FlipScvfMapType, numSD> scvfInfo_;
234};
235
236} // end namespace Dumux
237
238#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
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
the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
Definition: boundary/darcydarcy/couplingmapper.hh:49
bool isCoupledElement(Dune::index_constant< i >, std::size_t eIdx) const
Return if an element residual with index eIdx of domain i is coupled to domain j.
Definition: boundary/darcydarcy/couplingmapper.hh:190
void update(const CouplingManager &couplingManager)
Main update routine.
Definition: boundary/darcydarcy/couplingmapper.hh:81
const std::vector< std::size_t > & couplingStencil(Dune::index_constant< i > domainI, const std::size_t eIdxI, Dune::index_constant< j > domainJ) const
returns an iteratable container of all indices of degrees of freedom of domain j that couple with / i...
Definition: boundary/darcydarcy/couplingmapper.hh:175
std::size_t outsideElementIndex(Dune::index_constant< i > domainI, const SubControlVolumeFace< i > &scvf) const
Return the outside element index (the element index of the other domain)
Definition: boundary/darcydarcy/couplingmapper.hh:223
bool isCoupled(Dune::index_constant< i > domainI, const SubControlVolumeFace< i > &scvf) const
If the boundary entity is on a coupling boundary.
Definition: boundary/darcydarcy/couplingmapper.hh:199
std::size_t flipScvfIndex(Dune::index_constant< i > domainI, const SubControlVolumeFace< i > &scvf) const
Return the scvf index of the flipped scvf in the other domain.
Definition: boundary/darcydarcy/couplingmapper.hh:211
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
Algorithms that finds which geometric entites intersect.