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 |