GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/embedded/integrationpointsource.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 30 49 61.2%
Functions: 9 16 56.2%
Branches: 44 124 35.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-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 An integration point source class,
11 * i.e. sources located at a single point in space
12 * associated with a quadrature point
13 */
14
15 #ifndef DUMUX_INTEGRATION_POINTSOURCE_HH
16 #define DUMUX_INTEGRATION_POINTSOURCE_HH
17
18 #include <type_traits>
19 #include <dune/common/reservedvector.hh>
20 #include <dumux/common/pointsource.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/discretization/method.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup EmbeddedCoupling
28 * \brief An integration point source class with an identifier to attach data
29 * and a quadrature weight and integration element
30 */
31 template<class GlobalPosition, class SourceValues, typename IdType = std::size_t>
32 class IntegrationPointSource : public IdPointSource<GlobalPosition, SourceValues, IdType>
33 {
34 using ParentType = IdPointSource<GlobalPosition, SourceValues, IdType>;
35 using Scalar = std::decay_t<decltype(std::declval<SourceValues>()[0])>;
36
37 public:
38 //! Constructor for integration point sources
39 IntegrationPointSource(GlobalPosition pos, SourceValues values, IdType id,
40 Scalar qpweight, Scalar integrationElement,
41 std::size_t elementIndex)
42 : ParentType(pos, values, id)
43 , qpweight_(qpweight)
44 , integrationElement_(integrationElement)
45 , elementIndex_(elementIndex)
46 {}
47
48 //! Constructor for integration point sources, when there is no
49 //! value known at the time of initialization
50 170938 IntegrationPointSource(GlobalPosition pos, IdType id,
51 Scalar qpweight, Scalar integrationElement,
52 std::size_t elementIndex)
53 : ParentType(pos, id)
54 , qpweight_(qpweight)
55 , integrationElement_(integrationElement)
56 341876 , elementIndex_(elementIndex)
57 {}
58
59 Scalar quadratureWeight() const
60 { return qpweight_; }
61
62 Scalar integrationElement() const
63 { return integrationElement_; }
64
65 void setQuadratureWeight(const Scalar qpWeight)
66 { return qpweight_ = qpWeight; }
67
68 void setIntegrationElement(const Scalar ie)
69 { integrationElement_ = ie; }
70
71 //! The index of the element this integration point source is associated with
72 std::size_t elementIndex() const
73 { return elementIndex_; }
74
75 //! Convenience = operator overload modifying only the values
76 IntegrationPointSource& operator= (const SourceValues& values)
77 {
78 7048064 ParentType::operator=(values);
79 return *this;
80 }
81
82 //! Convenience = operator overload modifying only the values
83 IntegrationPointSource& operator= (Scalar s)
84 {
85 70016446 ParentType::operator=(s);
86 return *this;
87 }
88
89 private:
90 Scalar qpweight_;
91 Scalar integrationElement_;
92 std::size_t elementIndex_;
93 };
94
95 /*!
96 * \ingroup EmbeddedCoupling
97 * \brief A helper class calculating a DOF-index to point source map
98 */
99 class IntegrationPointSourceHelper
100 {
101
102 public:
103 //! calculate a DOF index to point source map from given vector of point sources
104 template<class GridGeometry, class PointSource, class PointSourceMap>
105 98 static void computePointSourceMap(const GridGeometry& gridGeometry,
106 const std::vector<PointSource>& sources,
107 PointSourceMap& pointSourceMap,
108 const std::string& paramGroup = "")
109 {
110
4/4
✓ Branch 0 taken 170938 times.
✓ Branch 1 taken 49 times.
✓ Branch 2 taken 170938 times.
✓ Branch 3 taken 49 times.
346674 for (const auto& source : sources)
111 {
112 // get the index of the element in which the point source falls
113 341876 const auto eIdx = source.elementIndex();
114 if constexpr (GridGeometry::discMethod == DiscretizationMethods::box
115 || GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
116 {
117
1/4
✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
13964 auto fvGeometry = localView(gridGeometry);
118 // check in which subcontrolvolume(s) we are
119
4/8
✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4494 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4494 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2252 times.
✗ Branch 11 not taken.
32884 const auto element = gridGeometry.boundingBoxTree().entitySet().entity(eIdx);
120
1/2
✓ Branch 1 taken 2252 times.
✗ Branch 2 not taken.
9460 fvGeometry.bindElement(element);
121 9460 const auto globalPos = source.position();
122
123
4/6
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4724 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
9460 static const bool boxPointSourceLumping = getParamFromGroup<bool>(paramGroup, "PointSource.EnableBoxLumping", true);
124
1/2
✓ Branch 0 taken 4730 times.
✗ Branch 1 not taken.
9460 if (boxPointSourceLumping)
125 {
126 // loop over all sub control volumes and check if the point source is inside
127 8988 constexpr int dim = GridGeometry::GridView::dimension;
128
2/2
✓ Branch 0 taken 236 times.
✓ Branch 1 taken 236 times.
9932 Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
129
4/4
✓ Branch 0 taken 27416 times.
✓ Branch 1 taken 4730 times.
✓ Branch 2 taken 27416 times.
✓ Branch 3 taken 4730 times.
73752 for (const auto& scv : scvs(fvGeometry))
130
3/4
✓ Branch 1 taken 27416 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6979 times.
✓ Branch 5 taken 20437 times.
54832 if (intersectsPointGeometry(globalPos, fvGeometry.geometry(scv)))
131 27916 scvIndices.push_back(scv.indexInElement());
132
133 // for all scvs that where tested positive add the point sources
134 // to the element/scv to point source map
135
2/2
✓ Branch 0 taken 6979 times.
✓ Branch 1 taken 4730 times.
42338 for (auto scvIdx : scvIndices)
136 {
137 13958 const auto key = std::make_pair(eIdx, scvIdx);
138 13958 if (pointSourceMap.count(key))
139
2/4
✓ Branch 1 taken 2063 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2063 times.
✗ Branch 5 not taken.
4126 pointSourceMap.at(key).push_back(source);
140 else
141
6/16
✓ Branch 1 taken 4916 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4916 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4916 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4916 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 4916 times.
✓ Branch 14 taken 4916 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
39328 pointSourceMap.insert({key, {source}});
142 // split equally on the number of matched scvs
143
1/2
✓ Branch 1 taken 6979 times.
✗ Branch 2 not taken.
13958 auto& s = pointSourceMap.at(key).back();
144 27916 s.setEmbeddings(scvIndices.size()*s.embeddings());
145 }
146 }
147
148 // distribute the sources using the basis function weights
149 else
150 {
151 const auto& localBasis = fvGeometry.feLocalBasis();
152 const auto ipLocal = element.geometry().local(globalPos);
153 using Scalar = std::decay_t<decltype(source.values()[0])>;
154 std::vector<typename Dune::FieldVector<Scalar, 1>> shapeValues;
155 localBasis.evaluateFunction(ipLocal, shapeValues);
156 for (const auto& scv : scvs(fvGeometry))
157 {
158 const auto key = std::make_pair(eIdx, scv.indexInElement());
159 if (pointSourceMap.count(key))
160 pointSourceMap.at(key).push_back(source);
161 else
162 pointSourceMap.insert({key, {source}});
163
164 // adjust the integration element
165 auto& s = pointSourceMap.at(key).back();
166 s.setIntegrationElement(shapeValues[scv.indexInElement()]*s.integrationElement());
167 }
168 }
169 }
170 else
171 {
172 // add the pointsource to the DOF map
173 332416 const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
174 664832 if (pointSourceMap.count(key))
175
3/6
✓ Branch 1 taken 149017 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 149017 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 149017 times.
✗ Branch 8 not taken.
596068 pointSourceMap.at(key).push_back(source);
176 else
177
6/16
✓ Branch 1 taken 17191 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17191 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 17191 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 17191 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 17191 times.
✓ Branch 14 taken 17191 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
137528 pointSourceMap.insert({key, {source}});
178 }
179 }
180 98 }
181 };
182
183 } // end namespace Dumux
184
185 #endif
186