version 3.10-dev
pointsource.hh
Go to the documentation of this file.
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//
14#ifndef DUMUX_POINTSOURCE_HH
15#define DUMUX_POINTSOURCE_HH
16
17#include <functional>
18
19#include <dune/common/reservedvector.hh>
26
28
29namespace Dumux {
30
37template<class PositionType, class ValueType>
39{
40public:
42 using Scalar = std::decay_t<decltype(std::declval<ValueType>()[0])>;
44 using GlobalPosition = PositionType;
46 using Values = ValueType;
47
50 : values_(values), pos_(pos), embeddings_(1) {}
51
53 // value known at the time of initialization
55 : values_(0.0), pos_(pos), embeddings_(1) {}
56
59 {
60 values_ += s;
61 return *this;
62 }
63
66 {
67 values_ -= s;
68 return *this;
69 }
70
73 {
74 values_ *= s;
75 return *this;
76 }
77
80 {
81 values_ /= s;
82 return *this;
83 }
84
87 {
89 return *this;
90 }
91
94 {
95 values_ = s;
96 return *this;
97 }
98
101 { return values_; }
102
105 { return pos_; }
106
108 // to the local residual in the problem in scvPointSources
109 // to be overloaded by derived classes
110 template<class Problem, class FVElementGeometry, class ElementVolumeVariables>
111 void update(const Problem &problem,
112 const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity &element,
113 const FVElementGeometry &fvGeometry,
114 const ElementVolumeVariables &elemVolVars,
115 const typename FVElementGeometry::SubControlVolume &scv)
116 {}
117
119 void setEmbeddings(std::size_t embeddings)
120 {
121 embeddings_ = embeddings;
122 }
123
135 std::size_t embeddings() const
136 {
137 return embeddings_;
138 }
139
140protected:
142private:
143 GlobalPosition pos_;
144 std::size_t embeddings_;
145};
146
154template<class GlobalPosition, class SourceValues, class I>
155class IdPointSource : public PointSource<GlobalPosition, SourceValues>
156{
158 using Scalar = typename ParentType::Scalar;
159
160public:
162 using IdType = I;
163
166 : ParentType(pos, values), id_(id) {}
167
169 // value known at the time of initialization
171 : ParentType(pos, SourceValues(0.0)), id_(id) {}
172
174 IdType id() const
175 { return id_; }
176
178 IdPointSource& operator= (const SourceValues& values)
179 {
181 return *this;
182 }
183
186 {
188 return *this;
189 }
190
191private:
192 IdType id_;
193};
194
199template<class TypeTag>
200class SolDependentPointSource : public PointSource<Dune::FieldVector<typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::ctype,
201 GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld>,
202 Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>>
203{
207 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
208 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
209 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
210 using Element = typename GridView::template Codim<0>::Entity;
211
212 static const int dimworld = GridView::dimensionworld;
213 using GlobalPosition = typename Dune::FieldVector<typename GridView::ctype, dimworld>;
214 // returns the PointSource values as PrimaryVariables
215 using ValueFunction = typename std::function<SourceValues(const Problem &problem,
216 const Element &element,
217 const FVElementGeometry &fvGeometry,
218 const ElementVolumeVariables &elemVolVars,
219 const SubControlVolume &scv)>;
220
222 using Scalar = typename ParentType::Scalar;
223
224public:
226 // value known at the time of initialization
227 SolDependentPointSource(GlobalPosition pos,
228 ValueFunction valueFunction)
229 : ParentType(pos, SourceValues(0.0)), valueFunction_(valueFunction) {}
230
232 // to the local residual in the problem in scvPointSources
233 // to be overloaded by derived classes
234 void update(const Problem &problem,
235 const Element &element,
236 const FVElementGeometry &fvGeometry,
237 const ElementVolumeVariables &elemVolVars,
238 const SubControlVolume &scv)
239 { this->values_ = valueFunction_(problem, element, fvGeometry, elemVolVars, scv); }
240
243 {
245 return *this;
246 }
247
250 {
252 return *this;
253 }
254
255private:
256 ValueFunction valueFunction_;
257};
258
266{
267public:
269 template<class GridGeometry, class PointSource, class PointSourceMap>
270 static void computePointSourceMap(const GridGeometry& gridGeometry,
271 const std::vector<PointSource>& sources,
272 PointSourceMap& pointSourceMap,
273 const std::string& paramGroup = "")
274 {
275 const auto& boundingBoxTree = gridGeometry.boundingBoxTree();
276
277 for (const auto& s : sources)
278 {
279 // compute in which elements the point source falls
280 const auto entities = intersectingEntities(s.position(), boundingBoxTree);
281
282 // continue with next point source if no intersection with the grid are found
283 if (entities.empty())
284 continue;
285
286 // make local copy of point source for the map
287 auto source = s;
288
289 // split the source values equally among all concerned entities
290 source.setEmbeddings(entities.size()*source.embeddings());
291
292 if constexpr (GridGeometry::discMethod == DiscretizationMethods::box
293 || GridGeometry::discMethod == DiscretizationMethods::fcdiamond)
294 {
295 // loop over all concerned elements
296 auto fvGeometry = localView(gridGeometry);
297 for (const auto eIdx : entities)
298 {
299 // check in which subcontrolvolume(s) we are
300 const auto element = boundingBoxTree.entitySet().entity(eIdx);
301 fvGeometry.bindElement(element);
302
303 const auto globalPos = source.position();
304 // loop over all sub control volumes and check if the point source is inside
305 constexpr int dim = GridGeometry::GridView::dimension;
306 Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
307 for (const auto& scv : scvs(fvGeometry))
308 if (intersectsPointGeometry(globalPos, fvGeometry.geometry(scv)))
309 scvIndices.push_back(scv.indexInElement());
310
311 // for all scvs that tested positive add the point sources
312 // to the element/scv to point source map
313 for (const auto scvIdx : scvIndices)
314 {
315 const auto key = std::make_pair(eIdx, scvIdx);
316 if (pointSourceMap.count(key))
317 pointSourceMap.at(key).push_back(source);
318 else
319 pointSourceMap.insert({key, {source}});
320 // split equally on the number of matched scvs
321 auto& s = pointSourceMap.at(key).back();
322 s.setEmbeddings(scvIndices.size()*s.embeddings());
323 }
324 }
325 }
326 else
327 {
328 for (const auto eIdx : entities)
329 {
330 // add the pointsource to the DOF map
331 const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0);
332 if (pointSourceMap.count(key))
333 pointSourceMap.at(key).push_back(source);
334 else
335 pointSourceMap.insert({key, {source}});
336 }
337 }
338 }
339 }
340};
341
342} // end namespace Dumux
343
344#endif
An axis-aligned bounding box volume hierarchy for dune grids.
A helper class calculating a sub control volume to point source map This class uses the bounding box ...
Definition: pointsource.hh:266
static void computePointSourceMap(const GridGeometry &gridGeometry, const std::vector< PointSource > &sources, PointSourceMap &pointSourceMap, const std::string &paramGroup="")
calculate a DOF index to point source map from given vector of point sources
Definition: pointsource.hh:270
A point source class with an identifier to attach data.
Definition: pointsource.hh:156
I IdType
export the id type
Definition: pointsource.hh:162
IdPointSource(GlobalPosition pos, IdType id)
Constructor for sol dependent point sources, when there is no.
Definition: pointsource.hh:170
IdPointSource(GlobalPosition pos, SourceValues values, IdType id)
Constructor for constant point sources.
Definition: pointsource.hh:165
IdType id() const
return the sources identifier
Definition: pointsource.hh:174
IdPointSource & operator=(const SourceValues &values)
Convenience = operator overload modifying only the values.
Definition: pointsource.hh:178
A point source base class.
Definition: pointsource.hh:39
void setEmbeddings(std::size_t embeddings)
set the number of embeddings for this point source
Definition: pointsource.hh:119
PointSource & operator*=(Scalar s)
Convenience *= operator overload modifying only the values.
Definition: pointsource.hh:72
ValueType Values
Export the value type.
Definition: pointsource.hh:46
PointSource & operator+=(Scalar s)
Convenience += operator overload modifying only the values.
Definition: pointsource.hh:58
PointSource(GlobalPosition pos)
Constructor for sol dependent point sources, when there is no.
Definition: pointsource.hh:54
std::decay_t< decltype(std::declval< ValueType >()[0])> Scalar
Export the scalar type.
Definition: pointsource.hh:42
Values values() const
return the source values
Definition: pointsource.hh:100
PositionType GlobalPosition
Export the position type.
Definition: pointsource.hh:44
Values values_
value of the point source for each equation
Definition: pointsource.hh:141
const GlobalPosition & position() const
return the source position
Definition: pointsource.hh:104
void update(const Problem &problem, const typename FVElementGeometry::GridGeometry::GridView::template Codim< 0 >::Entity &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const typename FVElementGeometry::SubControlVolume &scv)
an update function called before adding the value
Definition: pointsource.hh:111
PointSource & operator-=(Scalar s)
Convenience -= operator overload modifying only the values.
Definition: pointsource.hh:65
PointSource(GlobalPosition pos, Values values)
Constructor for constant point sources.
Definition: pointsource.hh:49
PointSource & operator=(const Values &values)
Convenience = operator overload modifying only the values.
Definition: pointsource.hh:86
PointSource & operator/=(Scalar s)
Convenience /= operator overload modifying only the values.
Definition: pointsource.hh:79
std::size_t embeddings() const
get the number of embeddings for this point source
Definition: pointsource.hh:135
A point source class for time dependent point sources.
Definition: pointsource.hh:203
SolDependentPointSource & operator=(const SourceValues &values)
Convenience = operator overload modifying only the values.
Definition: pointsource.hh:242
SolDependentPointSource(GlobalPosition pos, ValueFunction valueFunction)
Constructor for sol dependent point sources, when there is no.
Definition: pointsource.hh:227
void update(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolume &scv)
an update function called before adding the value
Definition: pointsource.hh:234
Defines all properties used in Dumux.
typename NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:34
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:26
bool intersectsPointGeometry(const Dune::FieldVector< ctype, dimworld > &point, const Geometry &g)
Find out whether a point is inside a three-dimensional geometry.
Definition: intersectspointgeometry.hh:28
std::vector< std::size_t > intersectingEntities(const Dune::FieldVector< ctype, dimworld > &point, const BoundingBoxTree< EntitySet > &tree, bool isCartesianGrid=false)
Compute all intersections between entities and a point.
Definition: intersectingentities.hh:102
typename GetProp< TypeTag, Property >::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:296
Algorithms that finds which geometric entities intersect.
Detect if a point intersects a geometry.
The available discretization methods in Dumux.
constexpr FCDiamond fcdiamond
Definition: method.hh:152
constexpr Box box
Definition: method.hh:147
Definition: adapt.hh:17
A helper to deduce a vector with the same size as numbers of equations.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.