3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
pointsourcedata.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
26#define DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH
27
28#include <vector>
29#include <dune/common/fvector.hh>
33
34namespace Dumux {
35
41template<class MDTraits>
43{
44 using Scalar = typename MDTraits::Scalar;
45 using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
46
47 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
48 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
49 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
50 template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
51 template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
52 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
53
54 static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
55 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
56
57 template<std::size_t id>
58 static constexpr bool isBox()
59 { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
60
61public:
62 void addBulkInterpolation(const ShapeValues& shapeValues,
63 const std::vector<GridIndex<bulkIdx>>& cornerIndices,
64 GridIndex<bulkIdx> eIdx)
65 {
66 static_assert(isBox<bulkIdx>(), "This interface is only available for the box method.");
67 bulkElementIdx_ = eIdx;
68 bulkCornerIndices_ = cornerIndices;
69 bulkShapeValues_ = shapeValues;
70 }
71
72 void addLowDimInterpolation(const ShapeValues& shapeValues,
73 const std::vector<GridIndex<lowDimIdx>>& cornerIndices,
74 GridIndex<lowDimIdx> eIdx)
75 {
76 static_assert(isBox<lowDimIdx>(), "This interface is only available for the box method.");
77 lowDimElementIdx_ = eIdx;
78 lowDimCornerIndices_ = cornerIndices;
79 lowDimShapeValues_ = shapeValues;
80 }
81
82 void addBulkInterpolation(GridIndex<bulkIdx> eIdx)
83 {
84 static_assert(!isBox<bulkIdx>(), "This interface is not available for the box method.");
85 bulkElementIdx_ = eIdx;
86 }
87
88 void addLowDimInterpolation(GridIndex<lowDimIdx> eIdx)
89 {
90 static_assert(!isBox<lowDimIdx>(), "This interface is not available for the box method.");
91 lowDimElementIdx_ = eIdx;
92 }
93
94 PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
95 {
96 PrimaryVariables<bulkIdx> bulkPriVars(0.0);
97 if (isBox<bulkIdx>())
98 {
99 for (int i = 0; i < bulkCornerIndices_.size(); ++i)
100 for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
101 bulkPriVars[priVarIdx] += sol[bulkCornerIndices_[i]][priVarIdx]*bulkShapeValues_[i];
102 }
103 else
104 {
105 bulkPriVars = sol[bulkElementIdx()];
106 }
107 return bulkPriVars;
108 }
109
110 PrimaryVariables<lowDimIdx> interpolateLowDim(const SolutionVector<lowDimIdx>& sol) const
111 {
112 PrimaryVariables<lowDimIdx> lowDimPriVars(0.0);
113 if (isBox<lowDimIdx>())
114 {
115 for (int i = 0; i < lowDimCornerIndices_.size(); ++i)
116 for (int priVarIdx = 0; priVarIdx < lowDimPriVars.size(); ++priVarIdx)
117 lowDimPriVars[priVarIdx] += sol[lowDimCornerIndices_[i]][priVarIdx]*lowDimShapeValues_[i];
118 }
119 else
120 {
121 lowDimPriVars = sol[lowDimElementIdx()];
122 }
123 return lowDimPriVars;
124 }
125
126 GridIndex<lowDimIdx> lowDimElementIdx() const
127 { return lowDimElementIdx_; }
128
129 GridIndex<bulkIdx> bulkElementIdx() const
130 { return bulkElementIdx_; }
131
132private:
133 ShapeValues bulkShapeValues_, lowDimShapeValues_;
134 std::vector<GridIndex<bulkIdx>> bulkCornerIndices_;
135 std::vector<GridIndex<lowDimIdx>> lowDimCornerIndices_;
136 GridIndex<bulkIdx> bulkElementIdx_;
137 GridIndex<lowDimIdx> lowDimElementIdx_;
138};
139
149template<class MDTraits>
151{
153 using Scalar = typename MDTraits::Scalar;
154 using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >;
155
156 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
157 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
158 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
159 template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex;
160 template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>;
161 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
162
163 static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
164 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
165
166 template<std::size_t id>
167 static constexpr bool isBox()
168 { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
169
170public:
171 PointSourceDataCircleAverage() : enableBulkCircleInterpolation_(false) {}
172
173 PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const
174 {
175 // bulk interpolation on the circle is only enabled for source in the
176 // lower dimensional domain if we use a circle distributed bulk sources
177 // (see coupling manager circle sources)
178 PrimaryVariables<bulkIdx> bulkPriVars(sol[0]);
179 bulkPriVars = 0.0;
180 if (enableBulkCircleInterpolation_)
181 {
182 // compute the average of the bulk privars over the circle around the integration point
183 // this computes $\bar{p} = \frac{1}{2\pi R} int_0^2*\pi p R \text{d}\theta.
184 if (isBox<bulkIdx>())
185 {
186 assert(circleCornerIndices_.size() == circleShapeValues_.size());
187
188 Scalar weightSum = 0.0;
189 for (std::size_t j = 0; j < circleStencil_.size(); ++j)
190 {
191 PrimaryVariables<bulkIdx> priVars(0.0);
192 const auto& cornerIndices = *(circleCornerIndices_[j]);
193 const auto& shapeValues = circleShapeValues_[j];
194 for (int i = 0; i < cornerIndices.size(); ++i)
195 {
196 const auto& localSol = sol[cornerIndices[i]];
197 const auto& shapeValue = shapeValues[i];
198 for (int priVarIdx = 0; priVarIdx < PrimaryVariables<bulkIdx>::size(); ++priVarIdx)
199 priVars[priVarIdx] += localSol[priVarIdx]*shapeValue;
200 }
201 // multiply with weight and add
202 priVars *= circleIpWeight_[j];
203 weightSum += circleIpWeight_[j];
204 bulkPriVars += priVars;
205 }
206 bulkPriVars /= weightSum;
207 }
208 else
209 {
210 Scalar weightSum = 0.0;
211 for (int j = 0; j < circleStencil_.size(); ++j)
212 {
213 for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx)
214 bulkPriVars[priVarIdx] += sol[circleStencil_[j]][priVarIdx]*circleIpWeight_[j];
215
216 weightSum += circleIpWeight_[j];
217 }
218 bulkPriVars /= weightSum;
219 }
220 return bulkPriVars;
221 }
222 else
223 {
224 return ParentType::interpolateBulk(sol);
225 }
226 }
227
228 void addCircleInterpolation(const std::vector<const std::vector<GridIndex<bulkIdx>>*>& circleCornerIndices,
229 const std::vector<ShapeValues>& circleShapeValues,
230 const std::vector<Scalar>& circleIpWeight,
231 const std::vector<GridIndex<bulkIdx>>& circleStencil)
232 {
233 circleCornerIndices_ = circleCornerIndices;
234 circleShapeValues_ = circleShapeValues;
235 circleIpWeight_ = circleIpWeight;
236 circleStencil_ = circleStencil;
237 enableBulkCircleInterpolation_ = true;
238
239 }
240
241 void addCircleInterpolation(const std::vector<Scalar>& circleIpWeight,
242 const std::vector<GridIndex<bulkIdx>>& circleStencil)
243 {
244 circleIpWeight_ = circleIpWeight;
245 circleStencil_ = circleStencil;
246 enableBulkCircleInterpolation_ = true;
247 }
248
249 const std::vector<GridIndex<bulkIdx>>& circleStencil() const
250 { return circleStencil_; }
251
252private:
253 std::vector<const std::vector<GridIndex<bulkIdx>>*> circleCornerIndices_;
254 std::vector<ShapeValues> circleShapeValues_;
255 std::vector<Scalar> circleIpWeight_;
256 std::vector<GridIndex<bulkIdx>> circleStencil_;
257 bool enableBulkCircleInterpolation_;
258};
259
260} // end namespace Dumux
261
262#endif
Defines the index types used for grid and local indices.
The available discretization methods in Dumux.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Struture to define the index types used for grid and local indices.
Definition: indextraits.hh:38
A vector of primary variables.
Definition: common/properties.hh:49
Vector containing all primary variable vector of the grid.
Definition: common/properties.hh:69
Definition: common/properties.hh:101
A point source data class used for integration in multidimension models.
Definition: pointsourcedata.hh:43
GridIndex< bulkIdx > bulkElementIdx() const
Definition: pointsourcedata.hh:129
void addLowDimInterpolation(const ShapeValues &shapeValues, const std::vector< GridIndex< lowDimIdx > > &cornerIndices, GridIndex< lowDimIdx > eIdx)
Definition: pointsourcedata.hh:72
GridIndex< lowDimIdx > lowDimElementIdx() const
Definition: pointsourcedata.hh:126
void addBulkInterpolation(GridIndex< bulkIdx > eIdx)
Definition: pointsourcedata.hh:82
PrimaryVariables< lowDimIdx > interpolateLowDim(const SolutionVector< lowDimIdx > &sol) const
Definition: pointsourcedata.hh:110
void addBulkInterpolation(const ShapeValues &shapeValues, const std::vector< GridIndex< bulkIdx > > &cornerIndices, GridIndex< bulkIdx > eIdx)
Definition: pointsourcedata.hh:62
PrimaryVariables< bulkIdx > interpolateBulk(const SolutionVector< bulkIdx > &sol) const
Definition: pointsourcedata.hh:94
void addLowDimInterpolation(GridIndex< lowDimIdx > eIdx)
Definition: pointsourcedata.hh:88
A point source data class used for integration in multidimension models.
Definition: pointsourcedata.hh:151
const std::vector< GridIndex< bulkIdx > > & circleStencil() const
Definition: pointsourcedata.hh:249
void addCircleInterpolation(const std::vector< Scalar > &circleIpWeight, const std::vector< GridIndex< bulkIdx > > &circleStencil)
Definition: pointsourcedata.hh:241
PointSourceDataCircleAverage()
Definition: pointsourcedata.hh:171
PrimaryVariables< bulkIdx > interpolateBulk(const SolutionVector< bulkIdx > &sol) const
Definition: pointsourcedata.hh:173
void addCircleInterpolation(const std::vector< const std::vector< GridIndex< bulkIdx > > * > &circleCornerIndices, const std::vector< ShapeValues > &circleShapeValues, const std::vector< Scalar > &circleIpWeight, const std::vector< GridIndex< bulkIdx > > &circleStencil)
Definition: pointsourcedata.hh:228
Declares all properties used in Dumux.