GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/discretization/projection/l2_projection.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 66 66 100.0%
Functions: 3 3 100.0%
Branches: 59 112 52.7%
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 L2-projections of analytic functions into a given function space
11 */
12 #ifndef DUMUX_DISCRETIZATION_L2_PROJECTION_HH
13 #define DUMUX_DISCRETIZATION_L2_PROJECTION_HH
14 #if HAVE_DUNE_FUNCTIONS
15
16 #include <vector>
17
18 #include <dune/common/fvector.hh>
19 #include <dune/common/fmatrix.hh>
20 #include <dune/common/parametertree.hh>
21 #include <dune/geometry/quadraturerules.hh>
22 #include <dune/istl/bcrsmatrix.hh>
23 #include <dune/istl/bvector.hh>
24 #include <dune/functions/gridfunctions/gridviewfunction.hh>
25
26 #include <dumux/linear/linearsolvertraits.hh>
27 #include <dumux/linear/linearalgebratraits.hh>
28 #include <dumux/linear/istlsolvers.hh>
29 #include <dumux/assembly/jacobianpattern.hh>
30
31 namespace Dumux {
32
33
34 template <class FEBasis>
35
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 class L2Projection
36 {
37 static constexpr int dim = FEBasis::GridView::dimension;
38 using FiniteElement = typename FEBasis::LocalView::Tree::FiniteElement;
39 using Scalar = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
40 using ShapeValue = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
41 using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
42 public:
43 using CoefficientVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
44
45 //! Parameters that can be passed to project()
46 struct Params
47 {
48 std::size_t maxIterations{100};
49 Scalar residualReduction{1e-13};
50 int verbosity{0};
51 };
52
53 1 L2Projection(const FEBasis& feBasis)
54 1 : feBasis_(feBasis)
55
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 , solver_()
56 {
57
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 3 solver_.setMatrix(std::make_shared<Matrix>(createMassMatrix_(feBasis)));
58 1 }
59
60 template <class Function>
61
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 CoefficientVector project(Function&& function, const Params& params = Params{}) const
62 {
63 1 CoefficientVector projection, rhs;
64
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 projection.resize(feBasis_.size());
65
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 rhs.resize(feBasis_.size());
66 1 rhs = 0.0;
67
68 // assemble right hand size
69 1 auto localFunc = localFunction(Dune::Functions::makeGridViewFunction(function, feBasis_.gridView()));
70
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto localView = feBasis_.localView();
71
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 101 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
 301 for (const auto& element : elements(feBasis_.gridView()))
72 {
73
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 localView.bind(element);
74 100 localFunc.bind(element);
75
76
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const auto& localFiniteElement = localView.tree().finiteElement();
77
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const int order = dim*localFiniteElement.localBasis().order();
78
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(element.type(), order);
79 100 const auto geometry = element.geometry();
80
81
2/2
✓ Branch 0 taken 900 times.
✓ Branch 1 taken 100 times.
 1000 for (auto&& qp : quad)
82 {
83 900 const auto weight = qp.weight();
84 1800 const auto ie = geometry.integrationElement(qp.position());
85 900 const auto globalPos = geometry.global(qp.position());
86
87 900 std::vector<ShapeValue> shapeValues;
88
1/4
✓ Branch 1 taken 900 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 1800 localFiniteElement.localBasis().evaluateFunction(geometry.local(globalPos), shapeValues);
89 900 const auto functionValue = localFunc(qp.position());
90
91
2/2
✓ Branch 1 taken 8100 times.
✓ Branch 2 taken 900 times.
 9000 for (int i = 0; i < localFiniteElement.localBasis().size(); ++i)
92 {
93 8100 const auto globalI = localView.index(i);
94 8100 rhs[globalI] += ie*weight*shapeValues[i]*functionValue;
95 }
96 }
97 }
98
99 // solve projection
100
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Dune::ParameterTree solverParams;
101
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 solverParams["maxit"] = std::to_string(params.maxIterations);
102
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 solverParams["reduction"] = std::to_string(params.residualReduction);
103
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 solverParams["verbose"] = std::to_string(params.verbosity);
104
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto solver = solver_; // copy the solver to modify the parameters
105
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 solver.setParams(solverParams);
106
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 solver.solve(projection, rhs);
107
108 1 return projection;
109 2 }
110
111 private:
112
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Matrix createMassMatrix_(const FEBasis& feBasis) const
113 {
114 1 Matrix massMatrix;
115
116
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto pattern = getFEJacobianPattern(feBasis);
117
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 pattern.exportIdx(massMatrix);
118
119
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto localView = feBasis.localView();
120
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 101 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 100 times.
✗ Branch 8 not taken.
 301 for (const auto& element : elements(feBasis.gridView()))
121 {
122
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 localView.bind(element);
123
124
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const auto& localFiniteElement = localView.tree().finiteElement();
125
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const int order = 2*(dim*localFiniteElement.localBasis().order()-1);
126
1/2
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
 100 const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(element.type(), order);
127 100 const auto geometry = element.geometry();
128
129
2/2
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 100 times.
 1700 for (auto&& qp : quad)
130 {
131 1600 const auto weight = qp.weight();
132 3200 const auto ie = geometry.integrationElement(qp.position());
133 1600 const auto globalPos = geometry.global(qp.position());
134
135 1600 std::vector<ShapeValue> shapeValues;
136
1/4
✓ Branch 1 taken 1600 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 3200 localFiniteElement.localBasis().evaluateFunction(geometry.local(globalPos), shapeValues);
137
138
2/2
✓ Branch 0 taken 14400 times.
✓ Branch 1 taken 1600 times.
 16000 for (int i = 0; i < localFiniteElement.localBasis().size(); ++i)
139 {
140
141
1/2
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
 14400 const auto globalI = localView.index(i);
142
1/2
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
 14400 massMatrix[globalI][globalI] += ie*weight*shapeValues[i]*shapeValues[i];
143
144
2/2
✓ Branch 0 taken 57600 times.
✓ Branch 1 taken 14400 times.
 72000 for (int j = i+1; j < localFiniteElement.localBasis().size(); ++j)
145 {
146
1/2
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
 57600 const auto globalJ = localView.index(j);
147
1/2
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
 57600 const auto value = ie*weight*shapeValues[i]*shapeValues[j];
148
2/4
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57600 times.
✗ Branch 5 not taken.
 57600 massMatrix[globalI][globalJ] += value;
149
1/2
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
 57600 massMatrix[globalJ][globalI] += value;
150 }
151 }
152 }
153 }
154
155 1 return massMatrix;
156 1 }
157
158 const FEBasis& feBasis_;
159 SSORCGIstlSolver<
160 SeqLinearSolverTraits, LinearAlgebraTraits<Matrix, CoefficientVector>
161 > solver_;
162 };
163
164 } // end namespace Dumux
165
166 #endif // HAVE_DUNE_FUNCTIONS
167 #endif
168