GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/discretization/projection/projector.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 195 229 85.2%
Functions: 15 15 100.0%
Branches: 215 518 41.5%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Discretization
10 * \brief Contains functionality for L2-projections from one function
11 * space into another, which can live both on the same or different
12 * grids of potentially different dimensionality.
13 * \note In the case of the function space bases living on the same grid
14 * an optimized implementation could be implemented. Currently, we
15 * require a glue object containing the intersections between two
16 * grids to be provided to the free factory functions.
17 */
18 #ifndef DUMUX_DISCRETIZATION_PROJECTOR_HH
19 #define DUMUX_DISCRETIZATION_PROJECTOR_HH
20
21 #include <algorithm>
22 #include <string>
23 #include <utility>
24 #include <type_traits>
25
26 #include <dune/common/timer.hh>
27 #include <dune/common/fmatrix.hh>
28 #include <dune/common/exceptions.hh>
29 #include <dune/common/promotiontraits.hh>
30 #include <dune/common/parametertree.hh>
31
32 #include <dune/geometry/quadraturerules.hh>
33 #include <dune/istl/matrixindexset.hh>
34 #include <dune/istl/bcrsmatrix.hh>
35 #include <dune/istl/bvector.hh>
36
37 #include <dumux/common/parameters.hh>
38 #include <dumux/linear/linearsolvertraits.hh>
39 #include <dumux/linear/linearalgebratraits.hh>
40 #include <dumux/linear/istlsolvers.hh>
41 #include <dumux/assembly/jacobianpattern.hh>
42 #include <dumux/discretization/functionspacebasis.hh>
43
44 namespace Dumux {
45
46 /*!
47 * \ingroup Discretization
48 * \brief Does an L2-projection from one discrete function space
49 * into another. The convenience functions makeProjectorPair
50 * or makeProjector can be used to create such a projection.
51 */
52 template<class ScalarType>
53 class Projector
54 {
55 using MatrixBlockType = Dune::FieldMatrix<ScalarType, 1, 1>;
56
57 public:
58 //! Export the scalar type
59 using Scalar = ScalarType;
60 //! Export the type of the projection matrices
61 using Matrix = Dune::BCRSMatrix< MatrixBlockType >;
62
63 //! Parameters that can be passed to project()
64 struct Params
65 {
66 std::size_t maxIterations{100};
67 Scalar residualReduction{1e-13};
68 int verbosity{0};
69 };
70
71 //! delete default constructor
72 Projector() = delete;
73
74 /*!
75 * \brief Constructor. Receives the mass and projection
76 * matrix that define the linear system describing
77 * the L2-projection from a function space into another.
78 */
79 29 Projector(Matrix&& massMatrix, Matrix&& projectionMatrix)
80
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 29 : massMat_(std::make_shared<Matrix>(std::move(massMatrix)))
81 29 , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix)))
82
3/6
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 29 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 29 times.
 58 , numDofsTarget_(massMat_->N())
83 {
84
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
 29 if (massMat_->N() != projMat_->N())
85 ✗ DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N());
86
87
1/2
✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
 29 massMatrixSolver_.setMatrix(massMat_);
88 29 }
89
90 /*!
91 * \brief Constructor for projection into a target space that occupies
92 * a larger geometric region than the domain space. In this case,
93 * the mass matrix can be chosen such that is is solved for only
94 * those dofs which will be populated with values from the projection.
95 * This requires an additional index map that maps the entries of the
96 * projected solution into the solution vector for the target space.
97 * Furthermore, the number of degrees of freedom must be specified to
98 * set up the coefficient vector with correct size for the target space.
99 */
100 1 Projector(Matrix&& massMatrix,
101 Matrix&& projectionMatrix,
102 std::vector<std::size_t>&& indexMap,
103 std::size_t numDofsTarget)
104
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 : massMat_(std::make_shared<Matrix>(std::move(massMatrix)))
105 1 , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix)))
106 1 , indexMapTarget_(std::move(indexMap))
107
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 2 , numDofsTarget_(numDofsTarget)
108 {
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (indexMapTarget_.size() != massMat_->N())
110 ✗ DUNE_THROW(Dune::InvalidStateException, "Target index map size mismatch: " << indexMapTarget_.size() << " vs " << massMat_->N());
111
112
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (massMat_->N() != projMat_->N())
113 ✗ DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N());
114
115
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (!indexMapTarget_.empty())
116
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (*std::max_element(indexMapTarget_.begin(), indexMapTarget_.end()) > numDofsTarget_)
117 ✗ DUNE_THROW(Dune::InvalidStateException, "Index map exceeds provided number of dofs in target domain!");
118
119
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 massMatrixSolver_.setMatrix(massMat_);
120 1 }
121
122 /*!
123 * \brief Project a solution u into up
124 * \param u The solution living on the domain space
125 * \param params Optional parameters for mass matrix solve
126 * \return The projection of u into the target space
127 */
128 template< class BlockType, std::enable_if_t<std::is_convertible<BlockType, ScalarType>::value, int> = 0 >
129 18 Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const
130 {
131 // be picky about size of u
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 if ( u.size() != projMat_->M())
133 ✗ DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch");
134
135 18 Dune::BlockVector<BlockType> up(massMat_->N());
136
137
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 auto rhs = up;
138 18 projMat_->mv(u, rhs);
139
140
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 Dune::ParameterTree solverParams;
141
3/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
 18 solverParams["maxit"] = std::to_string(params.maxIterations);
142
3/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
 18 solverParams["reduction"] = std::to_string(params.residualReduction);
143
2/4
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
 18 solverParams["verbose"] = std::to_string(params.verbosity);
144
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 auto solver = massMatrixSolver_; // copy the solver to modify the parameters
145
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 solver.setParams(solverParams);
146
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 solver.solve(up, rhs);
147
148 // if target space occupies a larger region, fill missing entries with zero
149
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 17 times.
 18 if (!indexMapTarget_.empty())
150 {
151
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Dune::BlockVector<BlockType> result(numDofsTarget_);
152
153 1 result = 0.0;
154
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 1 times.
 22 for (std::size_t i = 0; i < indexMapTarget_.size(); ++i)
155 21 result[indexMapTarget_[i]] = up[i];
156
157 1 return result;
158 1 }
159
160 17 return up;
161 25 }
162
163 /*!
164 * \brief Project a solution u into up
165 * \param u The solution living on the domain space
166 * \param params Optional parameters for mass matrix solve
167 * \return The projection of u into the target space
168 */
169 template< class BlockType, std::enable_if_t<!std::is_convertible<BlockType, ScalarType>::value, int> = 0 >
170 2 Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const
171 {
172 2 Dune::BlockVector<BlockType> result(numDofsTarget_);
173
174
3/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2 times.
 10 for (int pvIdx = 0; pvIdx < BlockType::size(); ++pvIdx)
175 {
176
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 Dune::BlockVector<Dune::FieldVector<Scalar, 1>> tmp(u.size());
177 11294 std::transform(u.begin(), u.end(), tmp.begin(), [pvIdx] (const auto& v) { return v[pvIdx]; });
178
179
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto p = project(tmp, params);
180
2/2
✓ Branch 0 taken 1124 times.
✓ Branch 1 taken 4 times.
 1128 for (std::size_t i = 0; i < p.size(); ++i)
181 1124 result[i][pvIdx] = p[i];
182 }
183
184 2 return result;
185 ✗ }
186
187 /*!
188 * \brief Returns the default parameters.
189 */
190 14 static Params defaultParams()
191
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 { return {}; }
192
193 private:
194 std::shared_ptr<Matrix> massMat_;
195 std::shared_ptr<Matrix> projMat_;
196
197 SSORCGIstlSolver<
198 SeqLinearSolverTraits, LinearAlgebraTraits<Matrix, Dune::BlockVector<Dune::FieldVector<Scalar, 1>>>
199 > massMatrixSolver_;
200
201 std::vector<std::size_t> indexMapTarget_;
202 std::size_t numDofsTarget_;
203 };
204
205 /*!
206 * \brief Traits class stating the type of projector between to bases
207 */
208 template<class FEBasisDomain, class FEBasisTarget>
209 class ProjectorTraits
210 {
211 using FiniteElementDomain = typename FEBasisDomain::LocalView::Tree::FiniteElement;
212 using FiniteElementTarget = typename FEBasisTarget::LocalView::Tree::FiniteElement;
213 using ScalarDomain = typename FiniteElementDomain::Traits::LocalBasisType::Traits::RangeFieldType;
214 using ScalarTarget = typename FiniteElementTarget::Traits::LocalBasisType::Traits::RangeFieldType;
215 using Scalar = typename Dune::PromotionTraits<ScalarDomain, ScalarTarget>::PromotedType;
216 public:
217 using Projector = Dumux::Projector< Scalar >;
218 };
219
220
221 // Projector construction details
222 namespace Detail {
223
224 /*!
225 * \brief Reduces a mass matrix and projection matrix such that they are composed
226 * of only those dofs that actually take part in the projection. Simultaneously,
227 * a container with the index map into the complete target space is filled so that
228 * the entries after projection can be assigned to the corresponding dof in the
229 * overall target space.
230 */
231 template<class Matrix>
232 1 void setupReducedMatrices(const Matrix& massMatrix, const Matrix& projMatrix, const std::vector<bool>& dofIsVoid,
233 Matrix& reducedM, Matrix& reducedP, std::vector<std::size_t>& expansionMap)
234 {
235 2 const std::size_t numNonVoidDofs = std::count_if(dofIsVoid.begin(), dofIsVoid.end(), [] (bool v) { return !v; });
236
237 // reduce matrices to only dofs that take part and create index map
238 1 std::vector<std::size_t> reductionMap(massMatrix.N());
239
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 expansionMap.resize(numNonVoidDofs);
240
241 std::size_t idxInReducedSpace = 0;
242
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
 442 for (std::size_t dofIdx = 0; dofIdx < dofIsVoid.size(); ++dofIdx)
243
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
 441 if (!dofIsVoid[dofIdx])
244 {
245 21 reductionMap[dofIdx] = idxInReducedSpace;
246 21 expansionMap[idxInReducedSpace] = dofIdx;
247 21 idxInReducedSpace++;
248 }
249
250 // create reduced M/P matrix
251 1 Dune::MatrixIndexSet patternMReduced, patternPReduced;
252
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 patternMReduced.resize(numNonVoidDofs, numNonVoidDofs);
253
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 patternPReduced.resize(numNonVoidDofs, projMatrix.M());
254
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
 442 for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
255
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 441 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
 441 if (!dofIsVoid[rowIt.index()])
256 {
257 21 const auto reducedRowIdx = reductionMap[rowIt.index()];
258
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 426 times.
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 21 times.
 426 for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
259
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 405 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 324 times.
 405 if (!dofIsVoid[colIt.index()])
260
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
 81 patternMReduced.add(reducedRowIdx, reductionMap[colIt.index()]);
261 }
262
263
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
 442 for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
264
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 441 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
 441 if (!dofIsVoid[rowIt.index()])
265 {
266 21 const auto reducedRowIdx = reductionMap[rowIt.index()];
267
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 156 times.
✓ Branch 2 taken 135 times.
✓ Branch 3 taken 21 times.
 156 for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
268
1/2
✓ Branch 1 taken 135 times.
✗ Branch 2 not taken.
 135 patternPReduced.add(reducedRowIdx, colIt.index());
269 }
270
271
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 patternMReduced.exportIdx(reducedM);
272
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 patternPReduced.exportIdx(reducedP);
273
274 // fill matrix entries
275
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
 442 for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
276
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 441 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
 441 if (!dofIsVoid[rowIt.index()])
277 {
278 21 const auto reducedRowIdx = reductionMap[rowIt.index()];
279
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 426 times.
✓ Branch 2 taken 405 times.
✓ Branch 3 taken 21 times.
 426 for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
280
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 405 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 324 times.
 405 if (!dofIsVoid[colIt.index()])
281
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
 81 reducedM[reducedRowIdx][reductionMap[colIt.index()]] = *colIt;
282 }
283
284
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
 442 for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
285
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 441 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
 441 if (!dofIsVoid[rowIt.index()])
286 {
287 21 const auto reducedRowIdx = reductionMap[rowIt.index()];
288
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 156 times.
✓ Branch 2 taken 135 times.
✓ Branch 3 taken 21 times.
 156 for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
289
1/2
✓ Branch 1 taken 135 times.
✗ Branch 2 not taken.
 135 reducedP[reducedRowIdx][colIt.index()] = *colIt;
290 }
291 1 }
292
293 /*!
294 * \brief Creates the matrices underlying l2-projections
295 * \tparam doBidirectional If false, the backward projection matrix is not assembled
296 * \param feBasisDomain The basis to the domain finite element space
297 * \param feBasisTarget The basis to the target finite element space
298 * \param glue The glue object containing the intersections between the two grids
299 * \param treatDiagonalZeroes If true, zero entries on the diagonal of the matrices
300 * that appear if the two domains occupy different geometric regions (and some
301 * dofs to not take part in the projection as a result) are substituted by ones.
302 * This substitution will lead to those dofs being mapped to zeroes in the target space.
303 * \returns An std::pair of projection matrices, where the first entry stores the
304 * matrices of the forward projection and the second entry stores those
305 * of the backward projection. The entries of the returned pair are itself
306 * std::pairs which store the mass matrix in the first and the projection
307 * matrix in the second entry.
308 */
309 template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType>
310 17 auto createProjectionMatrices(const FEBasisDomain& feBasisDomain,
311 const FEBasisTarget& feBasisTarget,
312 const GlueType& glue,
313 bool treatDiagonalZeroes = true)
314 {
315 // we assume that target dim <= domain dimension
316 static constexpr int domainDim = FEBasisDomain::GridView::dimension;
317 static constexpr int targetDim = FEBasisTarget::GridView::dimension;
318 static_assert(targetDim <= domainDim, "This expects target dim < domain dim, please swap arguments");
319
320 using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector;
321 using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector;
322
323 using ForwardProjectionMatrix = typename ForwardProjector::Matrix;
324 using BackwardProjectionMatrix = typename BackwardProjector::Matrix;
325
326
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 17 auto domainLocalView = feBasisDomain.localView();
327
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 17 auto targetLocalView = feBasisTarget.localView();
328
329 // determine mass matrix patterns (standard FE scheme pattern)
330 17 Dune::MatrixIndexSet backwardPatternM, forwardPatternM;
331
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 14 times.
✗ Branch 6 not taken.
 34 forwardPatternM = getFEJacobianPattern(feBasisTarget);
332
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 if (doBidirectional) backwardPatternM = getFEJacobianPattern(feBasisDomain);
333
334 // determine projection matrix patterns
335
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 17 Dune::MatrixIndexSet backwardPatternP, forwardPatternP;
336
3/5
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 17 forwardPatternP.resize(feBasisTarget.size(), feBasisDomain.size());
337
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 if (doBidirectional) backwardPatternP.resize(feBasisDomain.size(), feBasisTarget.size());
338
339 using std::max;
340 17 unsigned int maxBasisOrder = 0;
341
4/4
✓ Branch 1 taken 321334 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 321312 times.
✓ Branch 4 taken 14 times.
 335463 for (const auto& is : intersections(glue))
342 {
343 // since target dim <= domain dim there is maximum one!
344
1/2
✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
 335446 targetLocalView.bind( is.targetEntity(0) );
345 335446 const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
346
347
2/2
✓ Branch 0 taken 321356 times.
✓ Branch 1 taken 321334 times.
 670914 for (unsigned int nIdx = 0; nIdx < is.numDomainNeighbors(); ++nIdx)
348 {
349
1/2
✓ Branch 1 taken 321356 times.
✗ Branch 2 not taken.
 335468 domainLocalView.bind( is.domainEntity(nIdx) );
350
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 321356 times.
 335468 const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
351
352 // keep track of maximum basis order (used in integration)
353
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 321356 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 321347 times.
 335468 maxBasisOrder = max(maxBasisOrder, max(domainLocalBasis.order(), targetLocalBasis.order()));
354
355
2/2
✓ Branch 0 taken 860124 times.
✓ Branch 1 taken 321356 times.
 1287320 for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
356
2/2
✓ Branch 0 taken 3296820 times.
✓ Branch 1 taken 860124 times.
 4933104 for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
357 {
358
1/2
✓ Branch 1 taken 3296820 times.
✗ Branch 2 not taken.
 3981252 forwardPatternP.add(targetLocalView.index(j), domainLocalView.index(i));
359
1/2
✓ Branch 1 taken 1188 times.
✗ Branch 2 not taken.
 1188 if (doBidirectional) backwardPatternP.add(domainLocalView.index(i), targetLocalView.index(j));
360 }
361 }
362 }
363
364 // assemble matrices
365 17 ForwardProjectionMatrix forwardM, forwardP;
366
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 17 forwardPatternM.exportIdx(forwardM); forwardM = 0.0;
367
3/6
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 17 forwardPatternP.exportIdx(forwardP); forwardP = 0.0;
368
369 17 BackwardProjectionMatrix backwardM, backwardP;
370 if (doBidirectional)
371 {
372
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 backwardPatternM.exportIdx(backwardM); backwardM = 0.0;
373
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 backwardPatternP.exportIdx(backwardP); backwardP = 0.0;
374 }
375
376
4/4
✓ Branch 1 taken 321334 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 321312 times.
✓ Branch 4 taken 14 times.
 335463 for (const auto& is : intersections(glue))
377 {
378
1/2
✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
 335446 const auto& targetElement = is.targetEntity(0);
379 335446 const auto& targetElementGeometry = targetElement.geometry();
380
381
1/2
✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
 335446 targetLocalView.bind( targetElement );
382
1/2
✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
 335446 const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
383
384 // perform integration over intersection geometry
385 using IsGeometry = typename std::decay_t<decltype(is.geometry())>;
386 using ctype = typename IsGeometry::ctype;
387
388 335446 const auto& isGeometry = is.geometry();
389 335446 const int intOrder = maxBasisOrder + 1;
390
1/2
✓ Branch 2 taken 321334 times.
✗ Branch 3 not taken.
 335446 const auto& quad = Dune::QuadratureRules<ctype, IsGeometry::mydimension>::rule(isGeometry.type(), intOrder);
391
3/3
✓ Branch 0 taken 49436 times.
✓ Branch 1 taken 628534 times.
✓ Branch 2 taken 307200 times.
 1048674 for (auto&& qp : quad)
392 {
393 713228 const auto weight = qp.weight();
394 713228 const auto ie = isGeometry.integrationElement(qp.position());
395 713228 const auto globalPos = isGeometry.global(qp.position());
396
397
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 713228 std::vector< Dune::FieldVector<ctype, 1> > targetShapeVals;
398
1/4
✓ Branch 1 taken 663836 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 713228 targetLocalBasis.evaluateFunction(targetElementGeometry.local(globalPos), targetShapeVals);
399
400 // mass matrix entries target domain
401
2/2
✓ Branch 0 taken 2335620 times.
✓ Branch 1 taken 663836 times.
 3387536 for (unsigned int i = 0; i < targetLocalBasis.size(); ++i)
402 {
403
1/2
✓ Branch 1 taken 2335620 times.
✗ Branch 2 not taken.
 2674308 const auto dofIdxI = targetLocalView.index(i);
404
1/2
✓ Branch 1 taken 2335620 times.
✗ Branch 2 not taken.
 2674308 forwardM[dofIdxI][dofIdxI] += ie*weight*targetShapeVals[i]*targetShapeVals[i];
405
406
2/2
✓ Branch 0 taken 3908004 times.
✓ Branch 1 taken 2335620 times.
 7725384 for (unsigned int j = i+1; j < targetLocalBasis.size(); ++j)
407 {
408
1/2
✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
 5051076 const auto dofIdxJ = targetLocalView.index(j);
409
1/2
✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
 5051076 const auto value = ie*weight*targetShapeVals[i]*targetShapeVals[j];
410
2/4
✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3908004 times.
✗ Branch 5 not taken.
 5051076 forwardM[dofIdxI][dofIdxJ] += value;
411
1/2
✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
 5051076 forwardM[dofIdxJ][dofIdxI] += value;
412 }
413 }
414
415 // If targetDim < domainDim, there can be several "neighbors" if
416 // targetElement is aligned with a facet of domainElement. In
417 // this case make sure the basis functions are not added
418 // multiple times! (division by numNeighbors)
419 713228 const auto numNeighbors = is.numDomainNeighbors();
420
3/3
✓ Branch 0 taken 49480 times.
✓ Branch 1 taken 663836 times.
✓ Branch 2 taken 614400 times.
 1426500 for (unsigned int nIdx = 0; nIdx < numNeighbors; ++nIdx)
421 {
422 713272 const auto& domainElement = is.domainEntity(nIdx);
423
1/2
✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
 713272 domainLocalView.bind( domainElement );
424 713272 const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
425
426 713272 std::vector< Dune::FieldVector<ctype, 1> > domainShapeVals;
427
1/4
✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1426544 domainLocalBasis.evaluateFunction(domainElement.geometry().local(globalPos), domainShapeVals);
428
429 // add entries in matrices
430
2/2
✓ Branch 0 taken 2336280 times.
✓ Branch 1 taken 663880 times.
 3388240 for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
431 {
432
1/2
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
 2674968 const auto dofIdxDomain = domainLocalView.index(i);
433
1/2
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
 2674968 const auto domainShapeVal = domainShapeVals[i];
434 if (doBidirectional)
435 {
436
1/2
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
 792 backwardM[dofIdxDomain][dofIdxDomain] += ie*weight*domainShapeVal*domainShapeVal;
437
438
2/2
✓ Branch 0 taken 3168 times.
✓ Branch 1 taken 792 times.
 3960 for (unsigned int j = i+1; j < domainLocalBasis.size(); ++j)
439 {
440
1/2
✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
 3168 const auto dofIdxDomainJ = domainLocalView.index(j);
441
1/2
✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
 3168 const auto value = ie*weight*domainShapeVal*domainShapeVals[j];
442
2/4
✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3168 times.
✗ Branch 5 not taken.
 3168 backwardM[dofIdxDomain][dofIdxDomainJ] += value;
443
1/2
✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
 3168 backwardM[dofIdxDomainJ][dofIdxDomain] += value;
444 }
445 }
446
447
2/2
✓ Branch 0 taken 10153608 times.
✓ Branch 1 taken 2336280 times.
 15453408 for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
448 {
449
1/2
✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
 12778440 const auto dofIdxTarget = targetLocalView.index(j);
450
1/2
✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
 12778440 const auto entry = ie*weight*domainShapeVal*targetShapeVals[j];
451
452
2/4
✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2376 times.
✗ Branch 5 not taken.
 12778440 forwardP[dofIdxTarget][dofIdxDomain] += entry/numNeighbors;
453 if (doBidirectional)
454
1/2
✓ Branch 1 taken 2376 times.
✗ Branch 2 not taken.
 2376 backwardP[dofIdxDomain][dofIdxTarget] += entry;
455 }
456 }
457 }
458 }
459 }
460
461 // maybe treat zeroes on the diagonal
462
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 17 if (treatDiagonalZeroes)
463 {
464 ✗ for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget)
465 ✗ if (forwardM[dofIdxTarget][dofIdxTarget] == 0.0)
466 ✗ forwardM[dofIdxTarget][dofIdxTarget] = 1.0;
467
468 if (doBidirectional)
469 {
470 ✗ for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain)
471 ✗ if (backwardM[dofIdxDomain][dofIdxDomain] == 0.0)
472 ✗ backwardM[dofIdxDomain][dofIdxDomain] = 1.0;
473 }
474 }
475
476 return std::make_pair( std::make_pair(std::move(forwardM), std::move(forwardP)),
477
2/4
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
 51 std::make_pair(std::move(backwardM), std::move(backwardP)) );
478 17 }
479
480 /*!
481 * \brief Creates a projector class between two function space bases
482 * \tparam doBidirectional If false, the backward projection matrix is not assembled
483 * \returns an std::pair with the forward and backward projector
484 */
485 template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType>
486 17 auto makeProjectorPair(const FEBasisDomain& feBasisDomain,
487 const FEBasisTarget& feBasisTarget,
488 const GlueType& glue)
489 {
490 using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector;
491 using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector;
492
493 using ForwardProjectionMatrix = typename ForwardProjector::Matrix;
494 using BackwardProjectionMatrix = typename BackwardProjector::Matrix;
495
496 17 auto projectionMatrices = createProjectionMatrices<doBidirectional>(feBasisDomain, feBasisTarget, glue, false);
497 17 auto& forwardMatrices = projectionMatrices.first;
498 17 auto& backwardMatrices = projectionMatrices.second;
499
500 17 auto& forwardM = forwardMatrices.first;
501 17 auto& forwardP = forwardMatrices.second;
502
503 17 auto& backwardM = backwardMatrices.first;
504
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 17 auto& backwardP = backwardMatrices.second;
505
506 // determine the dofs that do not take part in intersections
507
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 17 std::vector<bool> isVoidTarget(forwardM.N(), false);
508
2/2
✓ Branch 0 taken 78355 times.
✓ Branch 1 taken 15 times.
 78934 for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget)
509
2/4
✓ Branch 1 taken 78355 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 78355 times.
 78917 if (forwardM[dofIdxTarget][dofIdxTarget] == 0.0)
510 ✗ isVoidTarget[dofIdxTarget] = true;
511
512
2/3
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 14 times.
 17 std::vector<bool> isVoidDomain;
513 if (doBidirectional)
514 {
515
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 isVoidDomain.resize(backwardM.N(), false);
516
2/2
✓ Branch 0 taken 441 times.
✓ Branch 1 taken 1 times.
 442 for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain)
517
3/4
✓ Branch 1 taken 441 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 420 times.
✓ Branch 4 taken 21 times.
 441 if (backwardM[dofIdxDomain][dofIdxDomain] == 0.0)
518 420 isVoidDomain[dofIdxDomain] = true;
519 }
520
521
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1 times.
 17 const bool hasVoidTarget = std::any_of(isVoidTarget.begin(), isVoidTarget.end(), [] (bool v) { return v; });
522 17 const bool hasVoidDomain = std::any_of(isVoidDomain.begin(), isVoidDomain.end(), [] (bool v) { return v; });
523
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1 times.
 17 if (!hasVoidDomain && !hasVoidTarget)
524 {
525
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 16 return std::make_pair(ForwardProjector(std::move(forwardM), std::move(forwardP)),
526
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 32 BackwardProjector(std::move(backwardM), std::move(backwardP)));
527 }
528
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 else if (!hasVoidDomain && hasVoidTarget)
529 {
530 ✗ std::vector<std::size_t> expansionMapTarget;
531 ✗ ForwardProjectionMatrix forwardMReduced, forwardPReduced;
532 ✗ setupReducedMatrices(forwardM, forwardP, isVoidTarget,
533 forwardMReduced, forwardPReduced, expansionMapTarget);
534
535 ✗ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
536 std::move(forwardPReduced),
537 std::move(expansionMapTarget),
538 forwardM.N()),
539 ✗ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
540 ✗ }
541
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 else if (hasVoidDomain && !hasVoidTarget)
542 {
543 if (doBidirectional)
544 {
545
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 std::vector<std::size_t> expansionMapDomain;
546 1 BackwardProjectionMatrix backwardMReduced, backwardPReduced;
547
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 setupReducedMatrices(backwardM, backwardP, isVoidDomain,
548 backwardMReduced, backwardPReduced, expansionMapDomain);
549
550
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
551
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 BackwardProjector(std::move(backwardMReduced),
552 std::move(backwardPReduced),
553 std::move(expansionMapDomain),
554 1 backwardM.N()) );
555 1 }
556 else
557 ✗ return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
558 ✗ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
559 }
560 else
561 {
562 ✗ std::vector<std::size_t> expansionMapTarget;
563 ✗ ForwardProjectionMatrix forwardMReduced, forwardPReduced;
564 ✗ setupReducedMatrices(forwardM, forwardP, isVoidTarget,
565 forwardMReduced, forwardPReduced, expansionMapTarget);
566
567 if (doBidirectional)
568 {
569 ✗ std::vector<std::size_t> expansionMapDomain;
570 ✗ BackwardProjectionMatrix backwardMReduced, backwardPReduced;
571 ✗ setupReducedMatrices(backwardM, backwardP, isVoidDomain,
572 backwardMReduced, backwardPReduced, expansionMapDomain);
573
574 ✗ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
575 std::move(forwardPReduced),
576 std::move(expansionMapTarget),
577 forwardM.N()),
578 ✗ BackwardProjector(std::move(backwardMReduced),
579 std::move(backwardPReduced),
580 std::move(expansionMapDomain),
581 ✗ backwardM.N()) );
582 ✗ }
583 else
584 ✗ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
585 std::move(forwardPReduced),
586 std::move(expansionMapTarget),
587 forwardM.N()),
588 ✗ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
589 ✗ }
590 17 }
591
592 } // end namespace Detail
593
594
595 /*!
596 * \ingroup Discretization
597 * \brief Creates a pair of projectors between the space with
598 * basis feBasisDomain to the space with basis feBasisTarget.
599 * \param feBasisDomain The domain finite element space basis
600 * \param feBasisTarget The target finite element space basis
601 * \param glue The glue object containing the intersections between the grids.
602 * \return An std::pair of projectors where the first is the forward
603 * projector from the space with basis feBasisDomain to
604 * the space with basis feBasisTarget and the second
605 * does the backward projection.
606 */
607 template< class FEBasisDomain, class FEBasisTarget, class GlueType >
608 1 auto makeProjectorPair(const FEBasisDomain& feBasisDomain,
609 const FEBasisTarget& feBasisTarget,
610 GlueType glue)
611 {
612 // we assume that target dim <= domain dimension
613 static constexpr int domainDim = FEBasisDomain::GridView::dimension;
614 static constexpr int targetDim = FEBasisTarget::GridView::dimension;
615 static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
616
617
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 return Detail::makeProjectorPair<true>(feBasisDomain, feBasisTarget, glue);
618 }
619
620 /*!
621 * \ingroup Discretization
622 * \brief Creates a forward projector from the space feBasisDomain
623 * to the space with basis feBasisTarget.
624 * \param feBasisDomain The domain finite element space basis
625 * \param feBasisTarget The target finite element space basis
626 * \param glue The glue object containing the intersections between the grids.
627 * \return The forward projector from the space with basis feBasisDomain
628 * to the space with basis feBasisTarget.
629 */
630 template< class FEBasisDomain, class FEBasisTarget, class GlueType >
631 14 auto makeProjector(const FEBasisDomain& feBasisDomain,
632 const FEBasisTarget& feBasisTarget,
633 GlueType glue)
634 {
635 // we assume that target dim <= domain dimension
636 static constexpr int domainDim = FEBasisDomain::GridView::dimension;
637 static constexpr int targetDim = FEBasisTarget::GridView::dimension;
638 static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
639
640 14 return Detail::makeProjectorPair<false>(feBasisDomain, feBasisTarget, glue).first;
641 }
642
643 /*!
644 * \brief Creates the matrices underlying l2-projections
645 * \param feBasisDomain The basis to the domain finite element space
646 * \param feBasisTarget The basis to the target finite element space
647 * \param glue The glue object containing the intersections between the two grids
648 * \returns An std::pair of projection matrices, where the first entry stores the
649 * matrices of the forward projection and the second entry stores those
650 * of the backward projection. The entries of the returned pair are itself
651 * std::pairs which store the mass matrix in the first and the projection
652 * matrix in the second entry.
653 */
654 template< class FEBasisDomain, class FEBasisTarget, class GlueType >
655 auto makeProjectionMatricesPair(const FEBasisDomain& feBasisDomain,
656 const FEBasisTarget& feBasisTarget,
657 GlueType glue)
658 {
659 // we assume that target dim <= domain dimension
660 static constexpr int domainDim = FEBasisDomain::GridView::dimension;
661 static constexpr int targetDim = FEBasisTarget::GridView::dimension;
662 static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments");
663
664 return Detail::createProjectionMatrices<true>(feBasisDomain, feBasisTarget, glue);
665 }
666
667 /*!
668 * \brief Creates the matrices underlying l2-projections
669 * \param feBasisDomain The basis to the domain finite element space
670 * \param feBasisTarget The basis to the target finite element space
671 * \param glue The glue object containing the intersections between the two grids
672 * \returns An std::pair of matrices, which store the mass matrix in the first and the
673 * projection matrix in the second entry.
674 */
675 template< class FEBasisDomain, class FEBasisTarget, class GlueType >
676 auto makeProjectionMatrices(const FEBasisDomain& feBasisDomain,
677 const FEBasisTarget& feBasisTarget,
678 GlueType glue)
679 {
680 // we assume that target dim <= domain dimension
681 static constexpr int domainDim = FEBasisDomain::GridView::dimension;
682 static constexpr int targetDim = FEBasisTarget::GridView::dimension;
683 static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments");
684
685 return Detail::createProjectionMatrices<false>(feBasisDomain, feBasisTarget, glue).first;
686 }
687
688 } // end namespace Dumux
689
690 #endif
691