GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 46 55 83.6%
Functions: 2 18 11.1%
Branches: 77 162 47.5%
Line Branch Exec Source
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 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup DarcyDarcyCoupling
10 * \copydoc Dumux::DarcyDarcyBoundaryCouplingMapper
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
14 #define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMAPPER_HH
15
16 #include <iostream>
17 #include <unordered_map>
18 #include <tuple>
19 #include <vector>
20
21 #include <dune/common/timer.hh>
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/indices.hh>
24
25 #include <dumux/geometry/intersectingentities.hh>
26 #include <dumux/discretization/method.hh>
27
28 namespace Dumux {
29
30 /*!
31 * \ingroup DarcyDarcyCoupling
32 * \brief the default mapper for conforming equal dimension boundary coupling between two domains (box or cc)
33 * \todo how to extend to arbitrary number of domains?
34 */
35 template<class MDTraits>
36 class DarcyDarcyBoundaryCouplingMapper
37 {
38 using Scalar = typename MDTraits::Scalar;
39
40 template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry;
41 template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace;
42 template<std::size_t i> using GridView = typename GridGeometry<i>::GridView;
43 template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity;
44
45 template<std::size_t i>
46 static constexpr auto domainIdx()
47 { return typename MDTraits::template SubDomain<i>::Index{}; }
48
49 template<std::size_t i>
50 static constexpr bool isCCTpfa()
51 { return GridGeometry<i>::discMethod == DiscretizationMethods::cctpfa; }
52
53 struct ScvfInfo
54 {
55 std::size_t eIdxOutside;
56 std::size_t flipScvfIdx;
57 };
58
59 using FlipScvfMapType = std::unordered_map<std::size_t, ScvfInfo>;
60 using MapType = std::unordered_map<std::size_t, std::vector<std::size_t>>;
61
62 static constexpr std::size_t numSD = MDTraits::numSubDomains;
63
64 public:
65 /*!
66 * \brief Main update routine
67 */
68 template<class CouplingManager>
69 17 void update(const CouplingManager& couplingManager)
70 {
71 // TODO: Box and multiple domains
72 static_assert(numSD == 2, "More than two subdomains not implemented!");
73 static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!");
74
75 17 Dune::Timer watch;
76 34 std::cout << "Initializing the coupling map..." << std::endl;
77
78
2/2
✓ Branch 0 taken 34 times.
✓ Branch 1 taken 17 times.
 51 for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
79 {
80 68 stencils_[domIdx].clear();
81 68 scvfInfo_[domIdx].clear();
82 }
83
84
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
 17 const auto& problem0 = couplingManager.problem(domainIdx<0>());
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
 17 const auto& problem1 = couplingManager.problem(domainIdx<1>());
86 17 const auto& gg0 = problem0.gridGeometry();
87 17 const auto& gg1 = problem1.gridGeometry();
88
89 48 isCoupledScvf_[0].resize(gg0.numScvf(), false);
90 34 isCoupledScvf_[1].resize(gg1.numScvf(), false);
91
92
2/4
✓ Branch 2 taken 136795 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 273590 for (const auto& element0 : elements(gg0.gridView()))
93 {
94
1/2
✓ Branch 1 taken 1578 times.
✗ Branch 2 not taken.
 138356 auto fvGeometry0 = localView(gg0);
95
1/2
✓ Branch 1 taken 1578 times.
✗ Branch 2 not taken.
 136778 fvGeometry0.bindElement(element0);
96
97
6/6
✓ Branch 0 taken 547112 times.
✓ Branch 1 taken 136778 times.
✓ Branch 2 taken 547112 times.
✓ Branch 3 taken 136778 times.
✓ Branch 4 taken 4640 times.
✓ Branch 5 taken 536160 times.
 820668 for (const auto& scvf0 : scvfs(fvGeometry0))
98 {
99 // skip all non-boundaries
100
2/2
✓ Branch 0 taken 4990 times.
✓ Branch 1 taken 542122 times.
 547112 if (!scvf0.boundary())
101 545822 continue;
102
103 // get elements intersecting with the scvf center
104 // for robustness add epsilon in unit outer normal direction
105 19960 const auto eps = (scvf0.ipGlobal() - element0.geometry().corner(0)).two_norm()*1e-8;
106 9980 auto globalPos = scvf0.ipGlobal(); globalPos.axpy(eps, scvf0.unitOuterNormal());
107
2/7
✓ Branch 1 taken 350 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 350 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 11270 const auto indices = intersectingEntities(globalPos, gg1.boundingBoxTree());
108
109 // skip if no intersection was found
110
4/4
✓ Branch 0 taken 3700 times.
✓ Branch 1 taken 1290 times.
✓ Branch 2 taken 3700 times.
✓ Branch 3 taken 1290 times.
 9980 if (indices.empty())
111
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3700 times.
 3700 continue;
112
113 // sanity check
114
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1290 times.
 2580 if (indices.size() > 1)
115 ✗ DUNE_THROW(Dune::InvalidStateException, "Are you sure your grids is conforming at the boundary?");
116
117 // add the pair to the multimap
118
2/4
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
 2580 const auto eIdx0 = gg0.elementMapper().index(element0);
119 1290 const auto eIdx1 = indices[0];
120
3/8
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1290 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2580 stencils_[0][eIdx0].push_back(eIdx1);
121
122 // mark the scvf and find and mark the flip scvf
123
3/6
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1290 times.
✗ Branch 8 not taken.
 3870 isCoupledScvf_[0][scvf0.index()] = true;
124
1/2
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
 1290 const auto& element1 = gg1.element(eIdx1);
125
2/4
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
 2580 auto fvGeometry1 = localView(gg1);
126
1/2
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
 1290 fvGeometry1.bindElement(element1);
127
128 using std::abs;
129
4/4
✓ Branch 0 taken 5160 times.
✓ Branch 1 taken 1290 times.
✓ Branch 2 taken 5160 times.
✓ Branch 3 taken 1290 times.
 7740 for (const auto& scvf1 : scvfs(fvGeometry1))
130
2/2
✓ Branch 0 taken 1376 times.
✓ Branch 1 taken 3784 times.
 5160 if (scvf1.boundary())
131
4/4
✓ Branch 0 taken 1290 times.
✓ Branch 1 taken 86 times.
✓ Branch 2 taken 1290 times.
✓ Branch 3 taken 86 times.
 2752 if (abs(scvf1.unitOuterNormal()*scvf0.unitOuterNormal() + 1) < eps)
132 {
133
3/6
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1290 times.
✗ Branch 8 not taken.
 3870 isCoupledScvf_[1][scvf1.index()] = true;
134
2/4
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
 2580 scvfInfo_[0][scvf0.index()] = ScvfInfo{eIdx1, scvf1.index()};
135
2/4
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
 2580 scvfInfo_[1][scvf1.index()] = ScvfInfo{eIdx0, scvf0.index()};
136 }
137 }
138 }
139
140 // create the inverse map for efficient access
141
2/2
✓ Branch 0 taken 1262 times.
✓ Branch 1 taken 17 times.
 1330 for (const auto& entry : stencils_[0])
142
4/4
✓ Branch 0 taken 1290 times.
✓ Branch 1 taken 1262 times.
✓ Branch 2 taken 1290 times.
✓ Branch 3 taken 1262 times.
 5076 for (const auto idx : entry.second)
143
3/6
✓ Branch 1 taken 1290 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1290 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1290 times.
✗ Branch 8 not taken.
 2580 stencils_[1][idx].push_back(entry.first);
144
145 51 std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
146 17 }
147
148 /*!
149 * \brief returns an iterable container of all indices of degrees of freedom of domain j
150 * that couple with / influence the element residual of the given element of domain i
151 *
152 * \param domainI the domain index of domain i
153 * \param eIdxI the index of the coupled element of domain í
154 * \param domainJ the domain index of domain j
155 *
156 * \note The element residual definition depends on the discretization scheme of domain i
157 * box: a container of the residuals of all sub control volumes
158 * cc : the residual of the (sub) control volume
159 * fem: the residual of the element
160 * \note This function has to be implemented by all coupling managers for all combinations of i and j
161 */
162 template<std::size_t i, std::size_t j>
163 ✗ const std::vector<std::size_t>& couplingStencil(Dune::index_constant<i> domainI,
164 const std::size_t eIdxI,
165 Dune::index_constant<j> domainJ) const
166 {
167 static_assert(i != j, "A domain cannot be coupled to itself!");
168 ✗ if (isCoupledElement(domainI, eIdxI))
169 ✗ return stencils_[i].at(eIdxI);
170 else
171 ✗ return emptyStencil_;
172 }
173
174 /*!
175 * \brief Return if an element residual with index eIdx of domain i is coupled to domain j
176 */
177 template<std::size_t i>
178 ✗ bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
179
4/10
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 6 taken 5440 times.
✓ Branch 7 taken 368040 times.
✓ Branch 10 taken 3160 times.
✓ Branch 11 taken 368040 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 1489360 { return static_cast<bool>(stencils_[i].count(eIdx)); }
180
181 /*!
182 * \brief If the boundary entity is on a coupling boundary
183 * \param domainI the domain index for which to compute the flux
184 * \param scvf the sub control volume face
185 */
186 template<std::size_t i>
187 ✗ bool isCoupled(Dune::index_constant<i> domainI,
188 const SubControlVolumeFace<i>& scvf) const
189 {
190
0/15
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 388652 return isCoupledScvf_[i].at(scvf.index());
191 }
192
193 /*!
194 * \brief Return the scvf index of the flipped scvf in the other domain
195 * \param domainI the domain index for which to compute the flux
196 * \param scvf the sub control volume face
197 */
198 template<std::size_t i>
199 ✗ std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
200 const SubControlVolumeFace<i>& scvf) const
201 {
202
8/14
✓ Branch 1 taken 54392 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54392 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5440 times.
✓ Branch 8 taken 51624 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3160 times.
✓ Branch 11 taken 51624 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 3160 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3160 times.
✗ Branch 17 not taken.
 218352 return scvfInfo_[i].at(scvf.index()).flipScvfIdx;
203 }
204
205 /*!
206 * \brief Return the outside element index (the element index of the other domain)
207 * \param domainI the domain index for which to compute the flux
208 * \param scvf the sub control volume face
209 */
210 template<std::size_t i>
211 ✗ std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
212 const SubControlVolumeFace<i>& scvf) const
213 {
214
2/4
✗ Branch 2 not taken.
✓ Branch 3 taken 5440 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3160 times.
 218352 return scvfInfo_[i].at(scvf.index()).eIdxOutside;
215 }
216
217 private:
218 std::array<MapType, numSD> stencils_;
219 std::vector<std::size_t> emptyStencil_;
220 std::array<std::vector<bool>, numSD> isCoupledScvf_;
221 std::array<FlipScvfMapType, numSD> scvfInfo_;
222 };
223
224 } // end namespace Dumux
225
226 #endif
227