GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/boundary/darcydarcy/couplingmapper.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 53 54 98.1%
Functions: 6 6 100.0%
Branches: 54 100 54.0%
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-FileCopyrightText: 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 18 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 18 Dune::Timer watch;
76 18 std::cout << "Initializing the coupling map..." << std::endl;
77
78
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 18 times.
 54 for (std::size_t domIdx = 0; domIdx < numSD; ++domIdx)
79 {
80 36 stencils_[domIdx].clear();
81 36 scvfInfo_[domIdx].clear();
82 }
83
84
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 const auto& problem0 = couplingManager.problem(domainIdx<0>());
85
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 const auto& problem1 = couplingManager.problem(domainIdx<1>());
86 18 const auto& gg0 = problem0.gridGeometry();
87 18 const auto& gg1 = problem1.gridGeometry();
88
89 18 isCoupledScvf_[0].resize(gg0.numScvf(), false);
90 18 isCoupledScvf_[1].resize(gg1.numScvf(), false);
91
92
3/5
✓ Branch 3 taken 1581 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1578 times.
✗ Branch 7 not taken.
✓ Branch 2 taken 135615 times.
 411552 for (const auto& element0 : elements(gg0.gridView()))
93 {
94
1/2
✓ Branch 1 taken 1578 times.
✗ Branch 2 not taken.
 137178 auto fvGeometry0 = localView(gg0);
95 137178 fvGeometry0.bindElement(element0);
96
97
4/4
✓ Branch 0 taken 543642 times.
✓ Branch 1 taken 5070 times.
✓ Branch 4 taken 548712 times.
✓ Branch 5 taken 137178 times.
 687200 for (const auto& scvf0 : scvfs(fvGeometry0))
98 {
99 // skip all non-boundaries
100
2/2
✓ Branch 0 taken 543642 times.
✓ Branch 1 taken 5070 times.
 548712 if (!scvf0.boundary())
101 547402 continue;
102
103 // get elements intersecting with the scvf center
104 // for robustness add epsilon in unit outer normal direction
105 15210 const auto eps = (scvf0.ipGlobal() - element0.geometry().corner(0)).two_norm()*1e-8;
106 5070 auto globalPos = scvf0.ipGlobal(); globalPos.axpy(eps, scvf0.unitOuterNormal());
107
1/2
✓ Branch 1 taken 5070 times.
✗ Branch 2 not taken.
 5070 const auto indices = intersectingEntities(globalPos, gg1.boundingBoxTree());
108
109 // skip if no intersection was found
110
2/2
✓ Branch 0 taken 3760 times.
✓ Branch 1 taken 1310 times.
 5070 if (indices.empty())
111 continue;
112
113 // sanity check
114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1310 times.
 1310 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.
 1310 const auto eIdx0 = gg0.elementMapper().index(element0);
119 1310 const auto eIdx1 = indices[0];
120
2/4
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1310 times.
✗ Branch 5 not taken.
 1310 stencils_[0][eIdx0].push_back(eIdx1);
121
122 // mark the scvf and find and mark the flip scvf
123
1/2
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
 1310 isCoupledScvf_[0][scvf0.index()] = true;
124
2/4
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1310 times.
✗ Branch 5 not taken.
 2620 const auto& element1 = gg1.element(eIdx1);
125
1/2
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
 1310 auto fvGeometry1 = localView(gg1);
126 1310 fvGeometry1.bindElement(element1);
127
128 using std::abs;
129
4/4
✓ Branch 0 taken 1398 times.
✓ Branch 1 taken 3842 times.
✓ Branch 2 taken 5240 times.
✓ Branch 3 taken 1310 times.
 6550 for (const auto& scvf1 : scvfs(fvGeometry1))
130
2/2
✓ Branch 0 taken 1398 times.
✓ Branch 1 taken 3842 times.
 5240 if (scvf1.boundary())
131
2/2
✓ Branch 0 taken 1310 times.
✓ Branch 1 taken 88 times.
 1398 if (abs(scvf1.unitOuterNormal()*scvf0.unitOuterNormal() + 1) < eps)
132 {
133
1/2
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
 1310 isCoupledScvf_[1][scvf1.index()] = true;
134
2/4
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1310 times.
✗ Branch 5 not taken.
 1310 scvfInfo_[0][scvf0.index()] = ScvfInfo{eIdx1, scvf1.index()};
135
1/2
✓ Branch 1 taken 1310 times.
✗ Branch 2 not taken.
 1310 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 1282 times.
✓ Branch 1 taken 18 times.
 1300 for (const auto& entry : stencils_[0])
142
2/2
✓ Branch 0 taken 1310 times.
✓ Branch 1 taken 1282 times.
 2592 for (const auto idx : entry.second)
143 1310 stencils_[1][idx].push_back(entry.first);
144
145 18 std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
146 18 }
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 2417616 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 4835232 if (isCoupledElement(domainI, eIdxI))
169 44152 return stencils_[i].at(eIdxI);
170 else
171 2373464 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 1984288 bool isCoupledElement(Dune::index_constant<i>, std::size_t eIdx) const
179
4/4
✓ Branch 1 taken 6600 times.
✓ Branch 2 taken 382480 times.
✓ Branch 4 taken 3920 times.
✓ Branch 5 taken 382480 times.
 3968576 { 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 824556 bool isCoupled(Dune::index_constant<i> domainI,
188 const SubControlVolumeFace<i>& scvf) const
189 {
190 824556 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 111096 std::size_t flipScvfIndex(Dune::index_constant<i> domainI,
200 const SubControlVolumeFace<i>& scvf) const
201 {
202
5/9
✓ Branch 1 taken 55552 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6600 times.
✓ Branch 5 taken 51624 times.
✓ Branch 7 taken 3920 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3920 times.
✗ Branch 11 not taken.
✗ Branch 6 not taken.
 111096 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 111096 std::size_t outsideElementIndex(Dune::index_constant<i> domainI,
212 const SubControlVolumeFace<i>& scvf) const
213 {
214
3/6
✓ Branch 1 taken 6600 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6600 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3920 times.
 111096 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