GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/embedded/couplingmanagerbase.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 176 195 90.3%
Functions: 272 593 45.9%
Branches: 441 779 56.6%

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 EmbeddedCoupling
10 * \brief Coupling manager for low-dimensional domains embedded in the bulk
11 * domain. Point sources on each integration point are computed by an AABB tree.
12 */
13
14 #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH
15 #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH
16
17 #include <iostream>
18 #include <fstream>
19 #include <string>
20 #include <utility>
21 #include <unordered_map>
22
23 #include <dune/common/timer.hh>
24 #include <dune/geometry/quadraturerules.hh>
25
26 #include <dumux/io/format.hh>
27 #include <dumux/parallel/parallel_for.hh>
28 #include <dumux/common/properties.hh>
29 #include <dumux/common/numeqvector.hh>
30 #include <dumux/assembly/coloring.hh>
31 #include <dumux/geometry/distance.hh>
32 #include <dumux/geometry/intersectingentities.hh>
33 #include <dumux/discretization/method.hh>
34 #include <dumux/multidomain/couplingmanager.hh>
35 #include <dumux/multidomain/glue.hh>
36 #include <dumux/multidomain/embedded/pointsourcedata.hh>
37 #include <dumux/multidomain/embedded/integrationpointsource.hh>
38 #include <dumux/multidomain/fvassembler.hh>
39
40 namespace Dumux {
41
42 //! the default point source traits
43 template<class MDTraits>
44 struct DefaultPointSourceTraits
45 {
46 private:
47 template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag;
48 template<std::size_t i> using GridGeometry = GetPropType<SubDomainTypeTag<i>, Properties::GridGeometry>;
49 template<std::size_t i> using NumEqVector = Dumux::NumEqVector<GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>>;
50 public:
51 //! export the point source type for domain i
52 template<std::size_t i>
53 using PointSource = IntegrationPointSource<typename GridGeometry<i>::GlobalCoordinate, NumEqVector<i>>;
54
55 //! export the point source helper type for domain i
56 template<std::size_t i>
57 using PointSourceHelper = IntegrationPointSourceHelper;
58
59 //! export the point source data type
60 using PointSourceData = Dumux::PointSourceData<MDTraits>;
61 };
62
63 /*!
64 * \ingroup EmbeddedCoupling
65 * \brief Manages the coupling between bulk elements and lower dimensional elements
66 * Point sources on each integration point are computed by an AABB tree.
67 */
68 template<class MDTraits, class Implementation, class PSTraits = DefaultPointSourceTraits<MDTraits>>
69 class EmbeddedCouplingManagerBase
70 : public CouplingManager<MDTraits>
71 {
72 using ParentType = CouplingManager<MDTraits>;
73 using Scalar = typename MDTraits::Scalar;
74 static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
75 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
76 using SolutionVector = typename MDTraits::SolutionVector;
77 using PointSourceData = typename PSTraits::PointSourceData;
78
79 // the sub domain type tags
80 template<std::size_t id> using PointSource = typename PSTraits::template PointSource<id>;
81 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
82 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
83 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
84 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
85 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
86 template<std::size_t id> using ElementMapper = typename GridGeometry<id>::ElementMapper;
87 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
88 template<std::size_t id> using ElementSeed = typename GridView<id>::Grid::template Codim<0>::EntitySeed;
89 template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
90 template<std::size_t id> using CouplingStencil = std::vector<GridIndex<id>>;
91
92 static constexpr int bulkDim = GridView<bulkIdx>::dimension;
93 static constexpr int lowDimDim = GridView<lowDimIdx>::dimension;
94 static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld;
95
96 template<std::size_t id>
97 static constexpr bool isBox()
98 { return GridGeometry<id>::discMethod == DiscretizationMethods::box; }
99
100 using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate;
101 using GlueType = MultiDomainGlue<GridView<bulkIdx>, GridView<lowDimIdx>, ElementMapper<bulkIdx>, ElementMapper<lowDimIdx>>;
102
103 public:
104 //! export traits
105 using MultiDomainTraits = MDTraits;
106 //! export the point source traits
107 using PointSourceTraits = PSTraits;
108 //! export stencil types
109 template<std::size_t id> using CouplingStencils = std::unordered_map<GridIndex<id>, CouplingStencil<id>>;
110
111 /*!
112 * \brief call this after grid adaption
113 */
114 void updateAfterGridAdaption(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry,
115 std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry)
116 {
117 glue_ = std::make_shared<GlueType>();
118 }
119
120 /*!
121 * \brief Constructor
122 */
123 27 EmbeddedCouplingManagerBase(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry,
124 std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry)
125 27 {
126
3/10
✓ Branch 3 taken 27 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 27 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 27 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
54 updateAfterGridAdaption(bulkGridGeometry, lowDimGridGeometry);
127 27 }
128
129 /*!
130 * \brief Methods to be accessed by main
131 */
132 // \{
133
134 27 void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem,
135 std::shared_ptr<Problem<lowDimIdx>> lowDimProblem,
136 const SolutionVector& curSol)
137 {
138 27 this->updateSolution(curSol);
139 54 this->setSubProblems(std::make_tuple(bulkProblem, lowDimProblem));
140
141 27 integrationOrder_ = getParam<int>("MixedDimension.IntegrationOrder", 1);
142 27 asImp_().computePointSourceData(integrationOrder_);
143 27 }
144
145 // \}
146
147 /*!
148 * \brief Methods to be accessed by the assembly
149 */
150 // \{
151
152 /*!
153 * \brief returns an iterable container of all indices of degrees of freedom of domain j
154 * that couple with / influence the element residual of the given element of domain i
155 *
156 * \param domainI the domain index of domain i
157 * \param element the coupled element of domain í
158 * \param domainJ the domain index of domain j
159 *
160 * \note The element residual definition depends on the discretization scheme of domain i
161 * box: a container of the residuals of all sub control volumes
162 * cc : the residual of the (sub) control volume
163 * fem: the residual of the element
164 * \note This function has to be implemented by all coupling managers for all combinations of i and j
165 */
166 template<std::size_t i, std::size_t j>
167 const CouplingStencil<j>& couplingStencil(Dune::index_constant<i> domainI,
168 const Element<i>& element,
169 Dune::index_constant<j> domainJ) const
170 {
171 static_assert(i != j, "A domain cannot be coupled to itself!");
172
173 const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
174 if (couplingStencils(domainI).count(eIdx))
175 return couplingStencils(domainI).at(eIdx);
176 else
177 return emptyStencil(domainI);
178 }
179
180 /*!
181 * \brief evaluates the element residual of a coupled element of domain i which depends on the variables
182 * at the degree of freedom with index dofIdxGlobalJ of domain j
183 *
184 * \param domainI the domain index of domain i
185 * \param localAssemblerI the local assembler assembling the element residual of an element of domain i
186 * \param domainJ the domain index of domain j
187 * \param dofIdxGlobalJ the index of the degree of freedom of domain j which has an influence on the element residual of domain i
188 *
189 * \note we only need to evaluate the source contribution to the residual here as the coupling term is the source
190 * \return the element residual
191 */
192 template<std::size_t i, std::size_t j, class LocalAssemblerI>
193 17727952 decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
194 const LocalAssemblerI& localAssemblerI,
195 Dune::index_constant<j> domainJ,
196 std::size_t dofIdxGlobalJ)
197 {
198 static_assert(i != j, "A domain cannot be coupled to itself!");
199
200 17727952 typename LocalAssemblerI::LocalResidual::ElementResidualVector residual;
201
202 17727952 const auto& element = localAssemblerI.element();
203 17727952 const auto& fvGeometry = localAssemblerI.fvGeometry();
204 17727952 const auto& curElemVolVars = localAssemblerI.curElemVolVars();
205
206 35455904 residual.resize(fvGeometry.numScv());
207
5/6
✓ Branch 0 taken 10211274 times.
✓ Branch 1 taken 8863976 times.
✓ Branch 2 taken 10211274 times.
✓ Branch 3 taken 8863976 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8566938 times.
76301000 for (const auto& scv : scvs(fvGeometry))
208 {
209
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10211274 times.
✓ Branch 3 taken 4916892 times.
✗ Branch 4 not taken.
20422548 auto couplingSource = this->problem(domainI).scvPointSources(element, fvGeometry, curElemVolVars, scv);
210
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10211274 times.
✓ Branch 3 taken 4916892 times.
✗ Branch 4 not taken.
21491304 couplingSource += this->problem(domainI).source(element, fvGeometry, curElemVolVars, scv);
211 44133768 couplingSource *= -GridGeometry<i>::Extrusion::volume(fvGeometry, scv)*curElemVolVars[scv].extrusionFactor();
212 40845096 residual[scv.indexInElement()] = couplingSource;
213 }
214 17727952 return residual;
215 }
216
217 // \}
218
219 /* \brief Compute integration point point sources and associated data
220 *
221 * This method uses grid glue to intersect the given grids. Over each intersection
222 * we later need to integrate a source term. This method places point sources
223 * at each quadrature point and provides the point source with the necessary
224 * information to compute integrals (quadrature weight and integration element)
225 * \param order The order of the quadrature rule for integration of sources over an intersection
226 * \param verbose If the point source computation is verbose
227 */
228 3 void computePointSourceData(std::size_t order = 1, bool verbose = false)
229 {
230 // initialize the maps
231 // do some logging and profiling
232 3 Dune::Timer watch;
233 6 std::cout << "Initializing the point sources..." << std::endl;
234
235 // clear all internal members like pointsource vectors and stencils
236 // initializes the point source id counter
237 3 clear();
238
239 // precompute the vertex indices for efficiency for the box method
240 3 this->precomputeVertexIndices(bulkIdx);
241 3 this->precomputeVertexIndices(lowDimIdx);
242
243
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
3 const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
244
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
3 const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
245
246 // intersect the bounding box trees
247 3 glueGrids();
248
249 3 pointSourceData_.reserve(this->glue().size());
250 3 averageDistanceToBulkCell_.reserve(this->glue().size());
251
4/4
✓ Branch 0 taken 22608 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 22608 times.
✓ Branch 3 taken 3 times.
22611 for (const auto& is : intersections(this->glue()))
252 {
253 // all inside elements are identical...
254 22608 const auto& inside = is.targetEntity(0);
255 // get the intersection geometry for integrating over it
256 22608 const auto intersectionGeometry = is.geometry();
257
258 // get the Gaussian quadrature rule for the local intersection
259 22608 const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order);
260 45216 const std::size_t lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside);
261
262 // iterate over all quadrature points
263
4/4
✓ Branch 0 taken 22628 times.
✓ Branch 1 taken 22608 times.
✓ Branch 2 taken 22628 times.
✓ Branch 3 taken 22608 times.
90452 for (auto&& qp : quad)
264 {
265 // compute the coupling stencils
266
4/4
✓ Branch 0 taken 22748 times.
✓ Branch 1 taken 22628 times.
✓ Branch 2 taken 22748 times.
✓ Branch 3 taken 22628 times.
68004 for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
267 {
268
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 const auto& outside = is.domainEntity(outsideIdx);
269 45496 const std::size_t bulkElementIdx = bulkGridGeometry.elementMapper().index(outside);
270
271 // each quadrature point will be a point source for the sub problem
272 45496 const auto globalPos = intersectionGeometry.global(qp.position());
273 22748 const auto id = idCounter_++;
274 22748 const auto qpweight = qp.weight();
275
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 const auto ie = intersectionGeometry.integrationElement(qp.position());
276 45496 pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, ie, bulkElementIdx);
277
2/4
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22748 times.
✗ Branch 5 not taken.
45496 pointSources(bulkIdx).back().setEmbeddings(is.numDomainNeighbors());
278 45496 pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, lowDimElementIdx);
279
2/4
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22748 times.
✗ Branch 5 not taken.
45496 pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors());
280
281 // pre compute additional data used for the evaluation of
282 // the actual solution dependent source term
283
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 PointSourceData psData;
284
285 if constexpr (isBox<lowDimIdx>())
286 {
287 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
288 const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry();
289 ShapeValues shapeValues;
290 this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues);
291 psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx);
292 }
293 else
294 {
295
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 psData.addLowDimInterpolation(lowDimElementIdx);
296 }
297
298 // add data needed to compute integral over the circle
299 if constexpr (isBox<bulkIdx>())
300 {
301 using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >;
302 const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry();
303 ShapeValues shapeValues;
304 this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues);
305 psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx);
306 }
307 else
308 {
309
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 psData.addBulkInterpolation(bulkElementIdx);
310 }
311
312 // publish point source data in the global vector
313
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 this->pointSourceData().emplace_back(std::move(psData));
314
315 // compute average distance to bulk cell
316
2/4
✓ Branch 3 taken 22748 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22748 times.
✗ Branch 7 not taken.
45496 averageDistanceToBulkCell_.push_back(averageDistancePointGeometry(globalPos, outside.geometry()));
317
318 // export the lowdim coupling stencil
319 // we insert all vertices / elements and make it unique later
320 if (isBox<bulkIdx>())
321 {
322 const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx);
323 this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(),
324 vertices.begin(), vertices.end());
325 }
326 else
327 {
328
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx);
329 }
330
331 // export the bulk coupling stencil
332 // we insert all vertices / elements and make it unique later
333 if (isBox<lowDimIdx>())
334 {
335 const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx);
336 this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(),
337 vertices.begin(), vertices.end());
338
339 }
340 else
341 {
342
1/2
✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
22748 this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx);
343 }
344 }
345 }
346 }
347
348 // make stencils unique
349 using namespace Dune::Hybrid;
350 27 forEach(integralRange(std::integral_constant<std::size_t, MDTraits::numSubDomains>{}), [&](const auto domainIdx)
351 {
352
4/4
✓ Branch 0 taken 616 times.
✓ Branch 1 taken 3 times.
✓ Branch 3 taken 3016 times.
✓ Branch 4 taken 3 times.
3644 for (auto&& stencil : this->couplingStencils(domainIdx))
353 {
354 10896 std::sort(stencil.second.begin(), stencil.second.end());
355 18160 stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end());
356 }
357 });
358
359 9 std::cout << "took " << watch.elapsed() << " seconds." << std::endl;
360 3 }
361
362 /*!
363 * \brief Methods to be accessed by the subproblems
364 */
365 // \{
366
367 //! Return a reference to the pointSource data
368 const PointSourceData& pointSourceData(std::size_t id) const
369
3/4
✓ Branch 0 taken 670548 times.
✓ Branch 1 taken 322387 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3175871 times.
332480975 { return pointSourceData_[id]; }
370
371 //! Return a reference to the bulk problem
372 template<std::size_t id>
373 const GridView<id>& gridView(Dune::index_constant<id> domainIdx) const
374
9/17
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 23 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 19 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✓ Branch 16 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 6 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 6 times.
112 { return this->problem(domainIdx).gridGeometry().gridView(); }
375
376 //! Return data for a bulk point source with the identifier id
377 PrimaryVariables<bulkIdx> bulkPriVars(std::size_t id) const
378
8/8
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1018714 times.
✓ Branch 2 taken 29 times.
✓ Branch 3 taken 1152049 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 464198 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 324923 times.
339074989 { return pointSourceData_[id].interpolateBulk(this->curSol(bulkIdx)); }
379
380 //! Return data for a low dim point source with the identifier id
381 PrimaryVariables<lowDimIdx> lowDimPriVars(std::size_t id) const
382
6/6
✓ Branch 0 taken 670572 times.
✓ Branch 1 taken 7375757 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 10769174 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 1367059 times.
183917453 { return pointSourceData_[id].interpolateLowDim(this->curSol(lowDimIdx)); }
383
384 //! return the average distance to the coupled bulk cell center
385 Scalar averageDistance(std::size_t id) const
386
8/8
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1012774 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 1012774 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 324923 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 324923 times.
2675452 { return averageDistanceToBulkCell_[id]; }
387
388 //! Return reference to bulk point sources
389 const std::vector<PointSource<bulkIdx>>& bulkPointSources() const
390
2/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
52 { return std::get<bulkIdx>(pointSources_); }
391
392 //! Return reference to low dim point sources
393 const std::vector<PointSource<lowDimIdx>>& lowDimPointSources() const
394
2/4
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
54 { return std::get<lowDimIdx>(pointSources_); }
395
396 //! Return the point source if domain i
397 template<std::size_t i>
398 const std::vector<PointSource<i>>& pointSources(Dune::index_constant<i> dom) const
399 { return std::get<i>(pointSources_); }
400
401 //! Return reference to bulk coupling stencil member of domain i
402 template<std::size_t i>
403 const CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom) const
404
0/8
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
2340 { return std::get<i>(couplingStencils_); }
405
406 //! Return reference to point source data vector member
407 const std::vector<PointSourceData>& pointSourceData() const
408
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 155548614 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 156846621 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1405896 times.
313801131 { return pointSourceData_; }
409
410 //! Return a reference to an empty stencil
411 template<std::size_t i>
412 const CouplingStencil<i>& emptyStencil(Dune::index_constant<i> dom) const
413
4/8
✓ Branch 0 taken 856990 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 856990 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 8424 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8424 times.
7144304 { return std::get<i>(emptyStencil_); }
414
415 /*!
416 * \brief the solution vector of the subproblem
417 * \param domainIdx The domain index
418 * \note in case of numeric differentiation the solution vector always carries the deflected solution
419 */
420 template<std::size_t i>
421 auto& curSol(Dune::index_constant<i> domainIdx)
422
3/6
✓ Branch 0 taken 808306 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 856990 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48684 times.
✗ Branch 5 not taken.
13157336 { return ParentType::curSol(domainIdx); }
423
424 /*!
425 * \brief the solution vector of the subproblem
426 * \param domainIdx The domain index
427 * \note in case of numeric differentiation the solution vector always carries the deflected solution
428 */
429 template<std::size_t i>
430 const auto& curSol(Dune::index_constant<i> domainIdx) const
431
16/16
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1018714 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 1018714 times.
✓ Branch 4 taken 24 times.
✓ Branch 5 taken 7375757 times.
✓ Branch 6 taken 24 times.
✓ Branch 7 taken 8603541 times.
✓ Branch 8 taken 10 times.
✓ Branch 9 taken 1691982 times.
✓ Branch 10 taken 10 times.
✓ Branch 11 taken 464198 times.
✓ Branch 12 taken 10 times.
✓ Branch 13 taken 464198 times.
✓ Branch 14 taken 10 times.
✓ Branch 15 taken 464198 times.
686412786 { return ParentType::curSol(domainIdx); }
432
433 /*!
434 * \brief Compute colors for multithreaded assembly
435 */
436 26 void computeColorsForAssembly()
437 {
438 26 Dune::Timer timer;
439
440 // start with elements connected within the subdomains
441
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
26 const auto& bulkGG = asImp_().problem(bulkIdx).gridGeometry();
442
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
26 const auto& lowDimGG = asImp_().problem(lowDimIdx).gridGeometry();
443
444 52 auto connectedElementsBulk = Detail::computeConnectedElements(bulkGG);
445
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
26 auto connectedElementsLowDim = Detail::computeConnectedElements(lowDimGG);
446
447 // to the elements from the own domain, we have to add the extended source stencil (non-local couplings)
448 // we consider this only in the bulk domain
449
14/18
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 953896 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 3198 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 3198 times.
✓ Branch 12 taken 8224 times.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 11422 times.
✓ Branch 15 taken 1 times.
✓ Branch 17 taken 8224 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 8224 times.
✗ Branch 21 not taken.
1912822 for (const auto& element : elements(bulkGG.gridView()))
450 {
451
2/4
✓ Branch 1 taken 11422 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11422 times.
✗ Branch 5 not taken.
1924192 const auto eIdx = bulkGG.elementMapper().index(element);
452 986320 const auto& extendedSourceStencil = asImp_().extendedSourceStencil(eIdx);
453
4/4
✓ Branch 0 taken 72633 times.
✓ Branch 1 taken 962096 times.
✓ Branch 2 taken 72633 times.
✓ Branch 3 taken 962096 times.
2934697 for (auto dofIdx : extendedSourceStencil)
454 {
455 72633 const auto& elems = connectedElementsBulk[dofIdx];
456
4/4
✓ Branch 4 taken 64364 times.
✓ Branch 5 taken 8269 times.
✓ Branch 6 taken 64364 times.
✓ Branch 7 taken 8269 times.
290532 if (std::find(elems.begin(), elems.end(), eIdx) == elems.end())
457
2/4
✓ Branch 1 taken 64364 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64364 times.
✗ Branch 5 not taken.
128728 connectedElementsBulk[dofIdx].push_back(eIdx);
458 }
459 }
460
461 // then we also need a similar container that tells us
462 // which elements _from the other domain_ modify our dofs.
463 // That, we can obtain via the coupling stencils
464 52 std::array<std::vector<std::vector<std::size_t>>, 2> connectedElementsCoupling;
465
466 using namespace Dune::Hybrid;
467
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
546 forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i)
468 {
469
4/16
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
104 connectedElementsCoupling[i()].resize(asImp_().problem(i).gridGeometry().numDofs());
470
0/4
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
4168792 forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto j)
471 {
472 if constexpr (i != j)
473 {
474
24/28
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✓ Branch 5 taken 953896 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3198 times.
✓ Branch 8 taken 950675 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 3222 times.
✓ Branch 11 taken 8224 times.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 8225 times.
✓ Branch 14 taken 6173 times.
✓ Branch 15 taken 24 times.
✓ Branch 16 taken 14396 times.
✓ Branch 17 taken 1762 times.
✓ Branch 18 taken 1 times.
✓ Branch 19 taken 16134 times.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 1738 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✓ Branch 31 taken 44 times.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 44 times.
✓ Branch 34 taken 1 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 44 times.
1937251 for (const auto& element : elements(asImp_().problem(j).gridGeometry().gridView()))
475 {
476
5/11
✗ Branch 0 not taken.
✓ Branch 1 taken 962096 times.
✓ Branch 3 taken 8224 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 16134 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8224 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 44 times.
970050 const auto eIdx = asImp_().problem(j).gridGeometry().elementMapper().index(element);
477
11/11
✓ Branch 1 taken 30559 times.
✓ Branch 2 taken 953872 times.
✓ Branch 3 taken 24503 times.
✓ Branch 4 taken 962096 times.
✓ Branch 5 taken 2168 times.
✓ Branch 6 taken 53343 times.
✓ Branch 7 taken 7910 times.
✓ Branch 8 taken 47580 times.
✓ Branch 9 taken 7954 times.
✓ Branch 10 taken 2461 times.
✓ Branch 11 taken 44 times.
1042133 for (auto dofIdx : asImp_().couplingStencil(j, element, i))
478 {
479 144166 const auto& elems = connectedElementsCoupling[i()][dofIdx];
480
4/8
✓ Branch 4 taken 24503 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24503 times.
✗ Branch 7 not taken.
✓ Branch 12 taken 47580 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 47580 times.
✗ Branch 15 not taken.
288332 if (std::find(elems.begin(), elems.end(), eIdx) == elems.end())
481
3/6
✓ Branch 1 taken 2168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2168 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2168 times.
✗ Branch 8 not taken.
216249 connectedElementsCoupling[i()][dofIdx].push_back(eIdx);
482 }
483 }
484 }
485 });
486 });
487
488
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
5804730 forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i)
489 {
490
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 26 times.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
52 const auto& gg = asImp_().problem(i).gridGeometry();
491
492
7/17
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 25 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 25 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 25 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 25 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 25 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
208 std::vector<int> colors(gg.gridView().size(0), -1);
493
494 // pre-reserve some memory for helper arrays to avoid reallocation
495
6/16
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 11 taken 26 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 26 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 26 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
156 std::vector<int> neighborColors; neighborColors.reserve(200);
496
7/13
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 25 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
156 std::vector<bool> colorUsed; colorUsed.reserve(200);
497
498
1/2
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
52 auto& elementSets = std::get<i>(elementSets_);
499
2/4
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
104 const auto& connectedElements = std::get<i>(std::tie(connectedElementsBulk, connectedElementsLowDim));
500
501
20/24
✓ Branch 2 taken 7954 times.
✓ Branch 3 taken 26 times.
✓ Branch 4 taken 7954 times.
✓ Branch 5 taken 26 times.
✓ Branch 6 taken 7928 times.
✓ Branch 7 taken 26 times.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 953896 times.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 3198 times.
✓ Branch 18 taken 950651 times.
✓ Branch 19 taken 3222 times.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 8224 times.
✓ Branch 22 taken 3199 times.
✓ Branch 23 taken 8224 times.
✓ Branch 24 taken 1 times.
✓ Branch 26 taken 8224 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 8224 times.
✗ Branch 30 not taken.
1937303 for (const auto& element : elements(gg.gridView()))
502 {
503 // compute neighbor colors based on discretization-dependent stencil
504
4/4
✓ Branch 0 taken 7928 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 962070 times.
✓ Branch 3 taken 26 times.
970050 neighborColors.clear();
505
2/4
✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
970050 Detail::addNeighborColors(gg, element, colors, connectedElements, neighborColors);
506
507 962096 if constexpr (i == bulkIdx)
508 {
509 // make sure we also add all elements in the extended source stencil
510 // Detail::addNeighborColors only adds direct neighbors
511 // add everyone that also modifies the extended source stencil dofs
512
2/4
✓ Branch 1 taken 11422 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11422 times.
✗ Branch 5 not taken.
1924192 const auto eIdx = bulkGG.elementMapper().index(element);
513 986320 const auto& extendedSourceStencil = asImp_().extendedSourceStencil(eIdx);
514
4/4
✓ Branch 0 taken 72633 times.
✓ Branch 1 taken 962096 times.
✓ Branch 2 taken 72633 times.
✓ Branch 3 taken 962096 times.
2934697 for (auto dofIdx : extendedSourceStencil)
515
4/4
✓ Branch 0 taken 2023554 times.
✓ Branch 1 taken 72633 times.
✓ Branch 2 taken 2023554 times.
✓ Branch 3 taken 72633 times.
2314086 for (const auto nIdx : connectedElementsBulk[dofIdx])
516
2/4
✓ Branch 1 taken 2023554 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2023554 times.
✗ Branch 5 not taken.
4047108 neighborColors.push_back(colors[nIdx]);
517 }
518
519 // add neighbor colors based on coupling stencil
520
2/4
✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
8059694 forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto j)
521 {
522 970050 if constexpr (i != j)
523 {
524 // add all elements that manipulate the same coupled dofs
525
8/8
✓ Branch 1 taken 47580 times.
✓ Branch 2 taken 7954 times.
✓ Branch 3 taken 47580 times.
✓ Branch 4 taken 7954 times.
✓ Branch 6 taken 24503 times.
✓ Branch 7 taken 962096 times.
✓ Branch 8 taken 24503 times.
✓ Branch 9 taken 962096 times.
1042133 for (auto dofIdx : asImp_().couplingStencil(i, element, j))
526
8/8
✓ Branch 0 taken 1062304 times.
✓ Branch 1 taken 47580 times.
✓ Branch 2 taken 1062304 times.
✓ Branch 3 taken 47580 times.
✓ Branch 4 taken 470335 times.
✓ Branch 5 taken 24503 times.
✓ Branch 6 taken 470335 times.
✓ Branch 7 taken 24503 times.
1893054 for (auto eIdx : connectedElementsCoupling[j][dofIdx])
527 3065278 neighborColors.push_back(colors[eIdx]);
528 }
529 });
530
531 // find smallest color (positive integer) not in neighborColors
532
2/4
✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
970050 const auto color = Detail::smallestAvailableColor(neighborColors, colorUsed);
533
534 // assign color to element
535
12/12
✓ Branch 0 taken 1720 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 7560 times.
✓ Branch 3 taken 394 times.
✓ Branch 4 taken 1720 times.
✓ Branch 5 taken 11437 times.
✓ Branch 6 taken 942116 times.
✓ Branch 7 taken 558 times.
✓ Branch 8 taken 19409 times.
✓ Branch 9 taken 13 times.
✓ Branch 10 taken 11258 times.
✓ Branch 11 taken 164 times.
1941835 colors[gg.elementMapper().index(element)] = color;
536
537 // add element to the set of elements with the same color
538
8/8
✓ Branch 0 taken 7560 times.
✓ Branch 1 taken 394 times.
✓ Branch 2 taken 7560 times.
✓ Branch 3 taken 394 times.
✓ Branch 4 taken 961361 times.
✓ Branch 5 taken 735 times.
✓ Branch 6 taken 961361 times.
✓ Branch 7 taken 735 times.
1940100 if (color < elementSets.size())
539
5/10
✓ Branch 2 taken 7560 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 7560 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 961361 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 8082 times.
✓ Branch 12 taken 953279 times.
✗ Branch 13 not taken.
977003 elementSets[color].push_back(element.seed());
540 else
541
14/28
✓ Branch 1 taken 394 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 394 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 394 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 394 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 394 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 735 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 142 times.
✓ Branch 21 taken 593 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 142 times.
✓ Branch 24 taken 593 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 142 times.
✓ Branch 27 taken 593 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 142 times.
✓ Branch 30 taken 593 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
3387 elementSets.push_back(std::vector<ElementSeed<i>>{ element.seed() });
542 }
543 });
544
545
3/8
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
26 std::cout << Fmt::format("Colored in {} seconds:\n", timer.elapsed());
546
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
364 forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i)
547 {
548
12/26
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 15 taken 25 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 25 times.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 25 times.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 25 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 25 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 25 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
208 std::cout << Fmt::format(
549 "-- {} elements in subdomain {} with {} colors\n",
550
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
52 asImp_().problem(i).gridGeometry().gridView().size(0),
551
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 26 times.
104 i(), std::get<i>(elementSets_).size()
552 );
553 });
554 26 }
555
556 /*!
557 * \brief Execute assembly kernel in parallel
558 *
559 * \param domainId the domain index of domain i
560 * \param assembleElement kernel function to execute for one element
561 */
562 template<std::size_t i, class AssembleElementFunc>
563 1168 void assembleMultithreaded(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
564 {
565
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 584 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 584 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 584 times.
3504 if (std::get<i>(elementSets_).empty())
566 DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!");
567
568 // make this element loop run in parallel
569 // for this we have to color the elements so that we don't get
570 // race conditions when writing into the global matrix or modifying grid variable caches
571 // each color can be assembled using multiple threads
572
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 584 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 584 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 584 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 584 times.
1168 const auto& grid = this->problem(domainId).gridGeometry().gridView().grid();
573
4/4
✓ Branch 0 taken 25556 times.
✓ Branch 1 taken 584 times.
✓ Branch 2 taken 25556 times.
✓ Branch 3 taken 584 times.
108064 for (const auto& elements : std::get<i>(elementSets_))
574 {
575
2/2
✓ Branch 2 taken 24657 times.
✓ Branch 3 taken 899 times.
5315500 Dumux::parallelFor(elements.size(), [&](const std::size_t n)
576 {
577
0/20
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
5871196 const auto element = grid.entity(elements[n]);
578
2/13
✓ Branch 1 taken 328960 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 328960 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
2606638 assembleElement(element);
579 });
580 }
581 1168 }
582
583 /*!
584 * \brief Extended source stencil (for the bulk domain)
585 * \note The default returns an empty stencil but some of the embedded
586 * coupling managers use this
587 */
588 const CouplingStencil<bulkIdx>& extendedSourceStencil(std::size_t eIdx) const
589 96896 { return asImp_().emptyStencil(bulkIdx); }
590
591 protected:
592 using ParentType::curSol;
593
594 //! computes the vertex indices per element for the box method
595 template<std::size_t id>
596 12 void precomputeVertexIndices(Dune::index_constant<id> domainIdx)
597 {
598 // fill helper structure for box discretization
599 if constexpr (isBox<domainIdx>())
600 {
601 12 this->vertexIndices(domainIdx).resize(gridView(domainIdx).size(0));
602
13/16
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 3 taken 15489 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 40 times.
✓ Branch 6 taken 15448 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 8224 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 8224 times.
✓ Branch 12 taken 1 times.
✓ Branch 14 taken 8224 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 8224 times.
✗ Branch 18 not taken.
111224 for (const auto& element : elements(gridView(domainIdx)))
603 {
604 47420 constexpr int dim = GridView<domainIdx>::dimension;
605
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 23710 times.
✓ Branch 3 taken 8224 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8224 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8224 times.
✗ Branch 10 not taken.
47420 const auto eIdx = this->problem(domainIdx).gridGeometry().elementMapper().index(element);
606
1/2
✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
47420 this->vertexIndices(domainIdx, eIdx).resize(element.subEntities(dim));
607
5/5
✓ Branch 0 taken 123648 times.
✓ Branch 1 taken 56606 times.
✓ Branch 2 taken 123648 times.
✓ Branch 3 taken 48382 times.
✓ Branch 4 taken 8224 times.
638776 for (int i = 0; i < element.subEntities(dim); ++i)
608
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 156544 times.
✓ Branch 3 taken 32896 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32896 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 32896 times.
✗ Branch 10 not taken.
313088 this->vertexIndices(domainIdx, eIdx)[i] = this->problem(domainIdx).gridGeometry().vertexMapper().subIndex(element, i, dim);
609 }
610 }
611 12 }
612
613 //! compute the shape function for a given point and geometry
614 template<std::size_t i, class FVGG, class Geometry, class ShapeValues>
615 81883 void getShapeValues(Dune::index_constant<i> domainI, const FVGG& gridGeometry, const Geometry& geo, const GlobalPosition& globalPos, ShapeValues& shapeValues)
616 {
617 if constexpr (FVGG::discMethod == DiscretizationMethods::box)
618 {
619 81883 const auto ipLocal = geo.local(globalPos);
620
1/12
✗ Branch 0 not taken.
✓ Branch 1 taken 62275 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
81883 const auto& localBasis = this->problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis();
621
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
81883 localBasis.evaluateFunction(ipLocal, shapeValues);
622 }
623 else
624 DUNE_THROW(Dune::InvalidStateException, "Shape values requested for other discretization than box!");
625 81883 }
626
627 //! Clear all internal data members
628 26 void clear()
629 {
630 78 pointSources(bulkIdx).clear();
631 78 pointSources(lowDimIdx).clear();
632 78 couplingStencils(bulkIdx).clear();
633 78 couplingStencils(lowDimIdx).clear();
634 78 vertexIndices(bulkIdx).clear();
635 78 vertexIndices(lowDimIdx).clear();
636
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
26 pointSourceData_.clear();
637
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
26 averageDistanceToBulkCell_.clear();
638
639 26 idCounter_ = 0;
640 26 }
641
642 //! compute the intersections between the two grids
643 26 void glueGrids()
644 {
645
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
26 const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry();
646
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
26 const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry();
647
648 // intersect the bounding box trees
649 104 glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree());
650 26 }
651
652 //! Return reference to point source data vector member
653 std::vector<PointSourceData>& pointSourceData()
654
4/7
✓ Branch 1 taken 61202 times.
✓ Branch 2 taken 1709 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
62919 { return pointSourceData_; }
655
656 //! Return reference to average distances to bulk cell
657 std::vector<Scalar>& averageDistanceToBulkCell()
658
2/3
✓ Branch 1 taken 10310 times.
✓ Branch 2 taken 380 times.
✗ Branch 3 not taken.
10695 { return averageDistanceToBulkCell_; }
659
660 //! Return the point source if domain i
661 template<std::size_t i>
662 std::vector<PointSource<i>>& pointSources(Dune::index_constant<i> dom)
663
21/32
✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 27 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 26 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 83407 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 83407 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 83407 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 83408 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 81704 times.
✓ Branch 22 taken 1 times.
✓ Branch 23 taken 6 times.
✓ Branch 24 taken 81698 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 48640 times.
✓ Branch 27 taken 33058 times.
✓ Branch 28 taken 48640 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 33058 times.
✗ Branch 31 not taken.
171388 { return std::get<i>(pointSources_); }
664
665 //! Return reference to bulk coupling stencil member of domain i
666 template<std::size_t i>
667 CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom)
668
10/18
✓ Branch 5 taken 90113 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 90113 times.
✓ Branch 9 taken 1709 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 90113 times.
✓ Branch 12 taken 1709 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 90113 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 58950 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 58950 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 236 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 236 times.
✗ Branch 27 not taken.
1010164 { return std::get<i>(couplingStencils_); }
669
670 //! Return a reference to the vertex indices
671 template<std::size_t i>
672 std::vector<GridIndex<i>>& vertexIndices(Dune::index_constant<i> dom, GridIndex<i> eIdx)
673
15/30
✓ Branch 4 taken 42753 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 62209 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 62209 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 30006 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10626 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 10626 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 236 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 236 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 76 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 76 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 76 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 76 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 76 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 76 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 76 times.
✗ Branch 47 not taken.
743128 { return std::get<i>(vertexIndices_)[eIdx]; }
674
675 //! Return a reference to the vertex indices container
676 template<std::size_t i>
677 std::vector<std::vector<GridIndex<i>>>& vertexIndices(Dune::index_constant<i> dom)
678
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
64 { return std::get<i>(vertexIndices_); }
679
680 const GlueType& glue() const
681 115 { return *glue_; }
682
683 //! Returns the implementation of the problem (i.e. static polymorphism)
684 Implementation &asImp_()
685 { return *static_cast<Implementation *>(this); }
686
687 //! \copydoc asImp_()
688 const Implementation &asImp_() const
689 { return *static_cast<const Implementation *>(this); }
690
691 //! id generator for point sources
692 std::size_t idCounter_ = 0;
693
694 private:
695
696 //! the point source in both domains
697 std::tuple<std::vector<PointSource<bulkIdx>>, std::vector<PointSource<lowDimIdx>>> pointSources_;
698 std::vector<PointSourceData> pointSourceData_;
699 std::vector<Scalar> averageDistanceToBulkCell_;
700
701 //! Stencil data
702 std::tuple<std::vector<std::vector<GridIndex<bulkIdx>>>,
703 std::vector<std::vector<GridIndex<lowDimIdx>>>> vertexIndices_;
704 std::tuple<CouplingStencils<bulkIdx>, CouplingStencils<lowDimIdx>> couplingStencils_;
705 std::tuple<CouplingStencil<bulkIdx>, CouplingStencil<lowDimIdx>> emptyStencil_;
706
707 //! The glue object
708 std::shared_ptr<GlueType> glue_;
709
710 //! integration order for coupling source
711 int integrationOrder_ = 1;
712
713 //! coloring for multithreaded assembly
714 std::tuple<
715 std::deque<std::vector<ElementSeed<bulkIdx>>>,
716 std::deque<std::vector<ElementSeed<lowDimIdx>>>
717 > elementSets_;
718 };
719
720 } // end namespace Dumux
721
722 #endif
723