GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/projection/l2_projection.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 31 63 49.2%
Functions: 2 3 66.7%
Branches: 56 236 23.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-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 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 : feBasis_(feBasis)
55
2/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1 , solver_()
56 {
57
4/10
✓ 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.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 1 solver_.setMatrix(std::make_shared<Matrix>(createMassMatrix_(feBasis)));
58 1 }
59
60 template <class Function>
61 ✗ CoefficientVector project(Function&& function, const Params& params = Params{}) const
62 {
63 ✗ CoefficientVector projection, rhs;
64 ✗ projection.resize(feBasis_.size());
65 ✗ rhs.resize(feBasis_.size());
66 ✗ rhs = 0.0;
67
68 // assemble right hand size
69 ✗ auto localFunc = localFunction(Dune::Functions::makeGridViewFunction(function, feBasis_.gridView()));
70 ✗ auto localView = feBasis_.localView();
71 ✗ for (const auto& element : elements(feBasis_.gridView()))
72 {
73 ✗ localView.bind(element);
74 ✗ localFunc.bind(element);
75
76 ✗ const auto& localFiniteElement = localView.tree().finiteElement();
77 ✗ const int order = dim*localFiniteElement.localBasis().order();
78 ✗ const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(element.type(), order);
79 ✗ const auto geometry = element.geometry();
80
81 ✗ for (auto&& qp : quad)
82 {
83 ✗ const auto weight = qp.weight();
84 ✗ const auto ie = geometry.integrationElement(qp.position());
85 ✗ const auto globalPos = geometry.global(qp.position());
86
87 ✗ std::vector<ShapeValue> shapeValues;
88 ✗ localFiniteElement.localBasis().evaluateFunction(geometry.local(globalPos), shapeValues);
89 ✗ const auto functionValue = localFunc(qp.position());
90
91 ✗ for (int i = 0; i < localFiniteElement.localBasis().size(); ++i)
92 {
93 ✗ const auto globalI = localView.index(i);
94 ✗ rhs[globalI] += ie*weight*shapeValues[i]*functionValue;
95 }
96 }
97 }
98
99 // solve projection
100 ✗ Dune::ParameterTree solverParams;
101 ✗ solverParams["maxit"] = std::to_string(params.maxIterations);
102 ✗ solverParams["reduction"] = std::to_string(params.residualReduction);
103 ✗ solverParams["verbose"] = std::to_string(params.verbosity);
104 ✗ auto solver = solver_; // copy the solver to modify the parameters
105 ✗ solver.setParams(solverParams);
106 ✗ solver.solve(projection, rhs);
107
108 ✗ return projection;
109 }
110
111 private:
112 1 Matrix createMassMatrix_(const FEBasis& feBasis) const
113 {
114
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 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.
 2 auto localView = feBasis.localView();
120
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 101 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 100 times.
✗ Branch 11 not taken.
 203 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 100 const auto& localFiniteElement = localView.tree().finiteElement();
125
2/4
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100 times.
✗ Branch 5 not taken.
 100 const int order = 2*(dim*localFiniteElement.localBasis().order()-1);
126 200 const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(element.type(), order);
127 100 const auto geometry = element.geometry();
128
129
4/4
✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 100 times.
✓ Branch 2 taken 1600 times.
✓ Branch 3 taken 100 times.
 3500 for (auto&& qp : quad)
130 {
131 1600 const auto weight = qp.weight();
132 3200 const auto ie = geometry.integrationElement(qp.position());
133 3200 const auto globalPos = geometry.global(qp.position());
134
135
2/4
✓ Branch 1 taken 1600 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1600 times.
✗ Branch 4 not taken.
 4800 std::vector<ShapeValue> shapeValues;
136
2/6
✓ Branch 1 taken 1600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1600 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3200 localFiniteElement.localBasis().evaluateFunction(geometry.local(globalPos), shapeValues);
137
138
4/6
✓ Branch 1 taken 16000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16000 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14400 times.
✓ Branch 7 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
7/14
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14400 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 14400 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 14400 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 14400 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 14400 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 14400 times.
✗ Branch 20 not taken.
 100800 massMatrix[globalI][globalI] += ie*weight*shapeValues[i]*shapeValues[i];
143
144
4/6
✓ Branch 1 taken 72000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72000 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14400 times.
✓ Branch 7 taken 57600 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
4/8
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 57600 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 57600 times.
✗ Branch 11 not taken.
 230400 const auto value = ie*weight*shapeValues[i]*shapeValues[j];
148
5/10
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 57600 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 57600 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 57600 times.
✗ Branch 14 not taken.
 230400 massMatrix[globalI][globalJ] += value;
149
4/8
✓ Branch 1 taken 57600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 57600 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 57600 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 57600 times.
✗ Branch 11 not taken.
 230400 massMatrix[globalJ][globalI] += value;
150 }
151 }
152 }
153 }
154
155 1 return massMatrix;
156 }
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