3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
fluxoversurface.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 *****************************************************************************/
24#ifndef DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
25#define DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
26
27#include <numeric>
28#include <functional>
29#include <type_traits>
30#include <vector>
31
32#include <dune/common/exceptions.hh>
33#include <dune/common/fvector.hh>
34#include <dune/geometry/type.hh>
35#include <dune/geometry/multilineargeometry.hh>
36#include <dune/geometry/referenceelements.hh>
37
41
42namespace Dumux {
43
48template<class GridVariables, class SolutionVector, class ModelTraits, class LocalResidual>
50{
51 using Scalar = typename GridVariables::Scalar;
52 using GridGeometry = typename GridVariables::GridGeometry;
53 using FVElementGeometry = typename GridGeometry::LocalView;
54 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
55 using GridView = typename GridGeometry::GridView;
56 using VolumeVariables = typename GridVariables::VolumeVariables;
57 using Element = typename GridView::template Codim<0>::Entity;
58
59 using CellCenterPrimaryVariables = std::decay_t<decltype(std::declval<SolutionVector>()[GridGeometry::cellCenterIdx()][0])>;
60
61 enum {
62 // Grid and world dimension
63 dim = GridView::dimension,
64 dimWorld = GridView::dimensionworld
65 };
66
67 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
68
69 static constexpr auto surfaceDim = dimWorld - 1;
70 using SurfaceGeometryType = Dune::MultiLinearGeometry< Scalar, surfaceDim, dimWorld >;
71
75 template<int mydim, int coordDim, class Scalar=double>
76 class SurfaceData
77 {
78 using Geo = SurfaceGeometryType;
79 public:
80
81 SurfaceData() {}
82
83 SurfaceData(Geo&& geo)
84 {
85 values_.resize(1);
86 geometries_.push_back(std::move(geo));
87 }
88
89 SurfaceData(std::vector<Geo>&& geo)
90 {
91 values_.resize(geo.size());
92 geometries_ = std::forward<decltype(geo)>(geo);
93 }
94
95 void addValue(int surfaceIdx, const CellCenterPrimaryVariables& value)
96 {
97 values_[surfaceIdx] += value;
98 }
99
100 void addSubSurface(Geo&& geo)
101 {
102 values_.emplace_back(0.0);
103 geometries_.push_back(std::move(geo));
104 }
105
106 auto& subSurfaces() const
107 {
108 return geometries_;
109 }
110
111 CellCenterPrimaryVariables& value(int surfaceIdx) const
112 {
113 return values_[surfaceIdx];
114 }
115
116 auto& values() const
117 {
118 return values_;
119 }
120
121 void printSurfaceBoundaries(int surfaceIdx) const
122 {
123 const auto& geometry = geometries_[surfaceIdx];
124 for(int i = 0; i < geometry.corners(); ++i)
125 std::cout << "(" << geometry.corner(i) << ") ";
126 }
127
128 void resetValues()
129 {
130 std::fill(values_.begin(), values_.end(), CellCenterPrimaryVariables(0.0));
131 }
132
133 private:
134
135 std::vector<Geo> geometries_;
136 std::vector<CellCenterPrimaryVariables> values_;
137 };
138
139public:
140
141 using SurfaceList = std::vector<SurfaceGeometryType>;
142
146 template<class Sol>
147 FluxOverSurface(const GridVariables& gridVariables,
148 const Sol& sol)
149 : gridVariables_(gridVariables),
150 sol_(sol)
151 {
152 static_assert(std::is_same<Sol, SolutionVector>::value, "Make sure that sol has the same type as SolutionVector."
153 "Use StaggeredVtkOutputModule<GridVariables, decltype(sol)> when calling the constructor.");
154
155 verbose_ = getParamFromGroup<bool>(problem_().paramGroup(), "FluxOverSurface.Verbose", false);
156 }
157
164 void addSurface(const std::string& name, SurfaceList&& surfaces)
165 {
166 surfaces_[name] = SurfaceData<surfaceDim,dim>(std::forward<decltype(surfaces)>(surfaces));
167 }
168
177 void addSurface(const std::string& name, const GlobalPosition& p0, const GlobalPosition& p1)
178 {
179 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1}));
180 }
181
192 void addSurface(const std::string& name,
193 const GlobalPosition& p0,
194 const GlobalPosition& p1,
195 const GlobalPosition& p2,
196 const GlobalPosition& p3)
197 {
198 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1, p2, p3}));
199 }
200
206 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 2>>& corners)
207 {
208 return SurfaceGeometryType(Dune::GeometryTypes::line, corners);
209 }
210
216 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 3>>& corners)
217 {
218 return makeDuneQuadrilaterial(corners);
219 }
220
225 {
226 auto fluxType = [this](const auto& element,
227 const auto& fvGeometry,
228 const auto& elemVolVars,
229 const auto& elemFaceVars,
230 const auto& scvf,
231 const auto& elemFluxVarsCache)
232 {
233 LocalResidual localResidual(&problem_());
234
235 if (scvf.boundary())
236 return localResidual.computeBoundaryFluxForCellCenter(problem_(), element, fvGeometry, scvf, elemVolVars, elemFaceVars, /*elemBcTypes*/{}, elemFluxVarsCache);
237 else
238 return localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
239 };
240
241 calculateFluxes(fluxType);
242 }
243
248 {
249 auto fluxType = [](const auto& element,
250 const auto& fvGeometry,
251 const auto& elemVolVars,
252 const auto& elemFaceVars,
253 const auto& scvf,
254 const auto& elemFluxVarsCache)
255 {
256 CellCenterPrimaryVariables result(0.0);
257 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
258 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
259 const Scalar extrusionFactor = harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
260 result[0] = elemFaceVars[scvf].velocitySelf() * scvf.area() * extrusionFactor * scvf.directionSign();
261 return result;
262 };
263
264 calculateFluxes(fluxType);
265 }
266
279 template<class FluxType>
280 void calculateFluxes(const FluxType& fluxType)
281 {
282 // make sure to reset all the values of the surfaces, in case this method has been called already before
283 for(auto&& surface : surfaces_)
284 surface.second.resetValues();
285
286 // make sure not to iterate over the same dofs twice
287 std::vector<bool> dofVisited(problem_().gridGeometry().numFaceDofs(), false);
288
289 auto elemVolVars = localView(gridVariables_.curGridVolVars());
290 auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache());
291 auto elemFaceVars = localView(gridVariables_.curGridFaceVars());
292
293 for(auto&& element : elements(problem_().gridGeometry().gridView()))
294 {
295 auto fvGeometry = localView(problem_().gridGeometry());
296 fvGeometry.bind(element);
297
298 elemVolVars.bind(element, fvGeometry, sol_);
299 elemFaceVars.bind(element, fvGeometry, sol_);
300 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
301
302 for(auto && scvf : scvfs(fvGeometry))
303 {
304 const auto dofIdx = scvf.dofIndex();
305 // do nothing of the dof was already visited
306 if(dofVisited[dofIdx])
307 continue;
308
309 dofVisited[dofIdx] = true;
310
311 // iterate through all surfaces and check if the flux at the given position
312 // should be accounted for in the respective surface
313 for(auto&& surface : surfaces_)
314 {
315 const auto& subSurfaces = surface.second.subSurfaces();
316
317 for(int surfaceIdx = 0; surfaceIdx < subSurfaces.size(); ++surfaceIdx)
318 {
319 if(intersectsPointGeometry(scvf.center(), subSurfaces[surfaceIdx]))
320 {
321 const auto result = fluxType(element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
322
323 surface.second.addValue(surfaceIdx, result);
324
325 if(verbose_)
326 {
327 std::cout << "Flux at face " << scvf.center() << " (" << surface.first << "): " << result;
328 std::cout << " (directionIndex: " << scvf.directionIndex() << "; surface boundaries: ";
329 surface.second.printSurfaceBoundaries(surfaceIdx);
330 std::cout << ", surfaceIdx " << surfaceIdx << ")" << std::endl;
331 }
332 }
333 }
334 }
335 }
336 }
337 }
338
344 auto& values(const std::string& name) const
345 {
346 return surfaces_.at(name).values();
347 }
348
354 auto netFlux(const std::string& name) const
355 {
356 const auto& surfaceResults = values(name);
357 return std::accumulate(surfaceResults.begin(), surfaceResults.end(), CellCenterPrimaryVariables(0.0));
358 }
359
360private:
361
362 const auto& problem_() const { return gridVariables_.curGridVolVars().problem(); }
363
364 std::map<std::string, SurfaceData<surfaceDim ,dim> > surfaces_;
365 const GridVariables& gridVariables_;
366 const SolutionVector& sol_;
367 bool verbose_;
368};
369
370} //end namespace
371
372#endif
Detect if a point intersects a geometry.
Create Dune geometries from user-specified points.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:68
auto makeDuneQuadrilaterial(const std::vector< Dune::FieldVector< CoordScalar, 3 > > &points)
Creates a dune quadrilateral geometry given 4 corner points.
Definition: makegeometry.hh:152
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:38
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Class used to calculate fluxes over surfaces. This only works for the staggered grid discretization.
Definition: fluxoversurface.hh:50
std::vector< SurfaceGeometryType > SurfaceList
Definition: fluxoversurface.hh:141
auto & values(const std::string &name) const
Return the fluxes of the individual sub surface of a given name.
Definition: fluxoversurface.hh:344
void addSurface(const std::string &name, const GlobalPosition &p0, const GlobalPosition &p1)
Add a surface under a given name, specifying the surface's corner points. This is a specialization fo...
Definition: fluxoversurface.hh:177
void calculateVolumeFluxes()
Calculate the volume fluxes over all surfaces.
Definition: fluxoversurface.hh:247
void addSurface(const std::string &name, SurfaceList &&surfaces)
Add a collection of sub surfaces under a given name.
Definition: fluxoversurface.hh:164
auto netFlux(const std::string &name) const
Return the cumulative net fluxes of a surface of a given name.
Definition: fluxoversurface.hh:354
static SurfaceGeometryType makeSurface(const std::vector< Dune::FieldVector< Scalar, 2 > > &corners)
Creates a geometrical surface object for (2D).
Definition: fluxoversurface.hh:206
void addSurface(const std::string &name, const GlobalPosition &p0, const GlobalPosition &p1, const GlobalPosition &p2, const GlobalPosition &p3)
Add a surface under a given name, specifying the surface's corner points. This is a specialization fo...
Definition: fluxoversurface.hh:192
void calculateMassOrMoleFluxes()
Calculate the mass or mole fluxes over all surfaces.
Definition: fluxoversurface.hh:224
FluxOverSurface(const GridVariables &gridVariables, const Sol &sol)
The constructor.
Definition: fluxoversurface.hh:147
static SurfaceGeometryType makeSurface(const std::vector< Dune::FieldVector< Scalar, 3 > > &corners)
Creates a geometrical surface object for (3D).
Definition: fluxoversurface.hh:216
void calculateFluxes(const FluxType &fluxType)
Calculate the fluxes over all surfaces for a given flux type.
Definition: fluxoversurface.hh:280