3.3.0
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
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 Extrusion = Extrusion_t<GridGeometry>;
56 using GridView = typename GridGeometry::GridView;
57 using VolumeVariables = typename GridVariables::VolumeVariables;
58 using Element = typename GridView::template Codim<0>::Entity;
59
60 using CellCenterPrimaryVariables = std::decay_t<decltype(std::declval<SolutionVector>()[GridGeometry::cellCenterIdx()][0])>;
61
62 enum {
63 // Grid and world dimension
64 dim = GridView::dimension,
65 dimWorld = GridView::dimensionworld
66 };
67
68 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
69
70 static constexpr auto surfaceDim = dimWorld - 1;
71 using SurfaceGeometryType = Dune::MultiLinearGeometry< Scalar, surfaceDim, dimWorld >;
72
76 template<int mydim, int coordDim, class Scalar=double>
77 class SurfaceData
78 {
79 using Geo = SurfaceGeometryType;
80 public:
81
82 SurfaceData() {}
83
84 SurfaceData(Geo&& geo)
85 {
86 values_.resize(1);
87 geometries_.push_back(std::move(geo));
88 }
89
90 SurfaceData(std::vector<Geo>&& geo)
91 {
92 values_.resize(geo.size());
93 geometries_ = std::forward<decltype(geo)>(geo);
94 }
95
96 void addValue(int surfaceIdx, const CellCenterPrimaryVariables& value)
97 {
98 values_[surfaceIdx] += value;
99 }
100
101 void addSubSurface(Geo&& geo)
102 {
103 values_.emplace_back(0.0);
104 geometries_.push_back(std::move(geo));
105 }
106
107 auto& subSurfaces() const
108 {
109 return geometries_;
110 }
111
112 CellCenterPrimaryVariables& value(int surfaceIdx) const
113 {
114 return values_[surfaceIdx];
115 }
116
117 auto& values() const
118 {
119 return values_;
120 }
121
122 void printSurfaceBoundaries(int surfaceIdx) const
123 {
124 const auto& geometry = geometries_[surfaceIdx];
125 for(int i = 0; i < geometry.corners(); ++i)
126 std::cout << "(" << geometry.corner(i) << ") ";
127 }
128
129 void resetValues()
130 {
131 std::fill(values_.begin(), values_.end(), CellCenterPrimaryVariables(0.0));
132 }
133
134 private:
135
136 std::vector<Geo> geometries_;
137 std::vector<CellCenterPrimaryVariables> values_;
138 };
139
140public:
141
142 using SurfaceList = std::vector<SurfaceGeometryType>;
143
147 template<class Sol>
148 FluxOverSurface(const GridVariables& gridVariables,
149 const Sol& sol)
150 : gridVariables_(gridVariables),
151 sol_(sol)
152 {
153 static_assert(std::is_same<Sol, SolutionVector>::value, "Make sure that sol has the same type as SolutionVector."
154 "Use StaggeredVtkOutputModule<GridVariables, decltype(sol)> when calling the constructor.");
155
156 verbose_ = getParamFromGroup<bool>(problem_().paramGroup(), "FluxOverSurface.Verbose", false);
157 }
158
165 void addSurface(const std::string& name, SurfaceList&& surfaces)
166 {
167 surfaces_[name] = SurfaceData<surfaceDim,dim>(std::forward<decltype(surfaces)>(surfaces));
168 }
169
178 void addSurface(const std::string& name, const GlobalPosition& p0, const GlobalPosition& p1)
179 {
180 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1}));
181 }
182
193 void addSurface(const std::string& name,
194 const GlobalPosition& p0,
195 const GlobalPosition& p1,
196 const GlobalPosition& p2,
197 const GlobalPosition& p3)
198 {
199 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1, p2, p3}));
200 }
201
207 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 2>>& corners)
208 {
209 return SurfaceGeometryType(Dune::GeometryTypes::line, corners);
210 }
211
217 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 3>>& corners)
218 {
219 return makeDuneQuadrilaterial(corners);
220 }
221
226 {
227 auto fluxType = [this](const auto& element,
228 const auto& fvGeometry,
229 const auto& elemVolVars,
230 const auto& elemFaceVars,
231 const auto& scvf,
232 const auto& elemFluxVarsCache)
233 {
234 LocalResidual localResidual(&problem_());
235
236 if (scvf.boundary())
237 return localResidual.computeBoundaryFluxForCellCenter(problem_(), element, fvGeometry, scvf, elemVolVars, elemFaceVars, /*elemBcTypes*/{}, elemFluxVarsCache);
238 else
239 return localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
240 };
241
242 calculateFluxes(fluxType);
243 }
244
249 {
250 auto fluxType = [](const auto& element,
251 const auto& fvGeometry,
252 const auto& elemVolVars,
253 const auto& elemFaceVars,
254 const auto& scvf,
255 const auto& elemFluxVarsCache)
256 {
257 CellCenterPrimaryVariables result(0.0);
258 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
259 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
260 const Scalar extrusionFactor = harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
261 result[0] = elemFaceVars[scvf].velocitySelf() * Extrusion::area(scvf) * extrusionFactor * scvf.directionSign();
262 return result;
263 };
264
265 calculateFluxes(fluxType);
266 }
267
280 template<class FluxType>
281 void calculateFluxes(const FluxType& fluxType)
282 {
283 // make sure to reset all the values of the surfaces, in case this method has been called already before
284 for(auto&& surface : surfaces_)
285 surface.second.resetValues();
286
287 // make sure not to iterate over the same dofs twice
288 std::vector<bool> dofVisited(problem_().gridGeometry().numFaceDofs(), false);
289
290 auto elemVolVars = localView(gridVariables_.curGridVolVars());
291 auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache());
292 auto elemFaceVars = localView(gridVariables_.curGridFaceVars());
293
294 for(auto&& element : elements(problem_().gridGeometry().gridView()))
295 {
296 auto fvGeometry = localView(problem_().gridGeometry());
297 fvGeometry.bind(element);
298
299 elemVolVars.bind(element, fvGeometry, sol_);
300 elemFaceVars.bind(element, fvGeometry, sol_);
301 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
302
303 for(auto && scvf : scvfs(fvGeometry))
304 {
305 const auto dofIdx = scvf.dofIndex();
306 // do nothing of the dof was already visited
307 if(dofVisited[dofIdx])
308 continue;
309
310 dofVisited[dofIdx] = true;
311
312 // iterate through all surfaces and check if the flux at the given position
313 // should be accounted for in the respective surface
314 for(auto&& surface : surfaces_)
315 {
316 const auto& subSurfaces = surface.second.subSurfaces();
317
318 for(int surfaceIdx = 0; surfaceIdx < subSurfaces.size(); ++surfaceIdx)
319 {
320 if(intersectsPointGeometry(scvf.center(), subSurfaces[surfaceIdx]))
321 {
322 const auto result = fluxType(element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
323
324 surface.second.addValue(surfaceIdx, result);
325
326 if(verbose_)
327 {
328 std::cout << "Flux at face " << scvf.center() << " (" << surface.first << "): " << result;
329 std::cout << " (directionIndex: " << scvf.directionIndex() << "; surface boundaries: ";
330 surface.second.printSurfaceBoundaries(surfaceIdx);
331 std::cout << ", surfaceIdx " << surfaceIdx << ")" << std::endl;
332 }
333 }
334 }
335 }
336 }
337 }
338 }
339
345 auto& values(const std::string& name) const
346 {
347 return surfaces_.at(name).values();
348 }
349
355 auto netFlux(const std::string& name) const
356 {
357 const auto& surfaceResults = values(name);
358 return std::accumulate(surfaceResults.begin(), surfaceResults.end(), CellCenterPrimaryVariables(0.0));
359 }
360
361private:
362
363 const auto& problem_() const { return gridVariables_.curGridVolVars().problem(); }
364
365 std::map<std::string, SurfaceData<surfaceDim ,dim> > surfaces_;
366 const GridVariables& gridVariables_;
367 const SolutionVector& sol_;
368 bool verbose_;
369};
370
371} //end namespace
372
373#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Helper classes to compute the integration elements.
auto makeDuneQuadrilaterial(const std::vector< Dune::FieldVector< CoordScalar, 3 > > &points)
Creates a dune quadrilateral geometry given 4 corner points.
Definition: geometry/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: geometry/intersectspointgeometry.hh:38
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
Definition: adapt.hh:29
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
constexpr Line line
Definition: couplingmanager1d3d_line.hh:52
constexpr Surface surface
Definition: couplingmanager1d3d_surface.hh:57
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:142
auto & values(const std::string &name) const
Return the fluxes of the individual sub surface of a given name.
Definition: fluxoversurface.hh:345
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:178
void calculateVolumeFluxes()
Calculate the volume fluxes over all surfaces.
Definition: fluxoversurface.hh:248
void addSurface(const std::string &name, SurfaceList &&surfaces)
Add a collection of sub surfaces under a given name.
Definition: fluxoversurface.hh:165
auto netFlux(const std::string &name) const
Return the cumulative net fluxes of a surface of a given name.
Definition: fluxoversurface.hh:355
static SurfaceGeometryType makeSurface(const std::vector< Dune::FieldVector< Scalar, 2 > > &corners)
Creates a geometrical surface object for (2D).
Definition: fluxoversurface.hh:207
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:193
void calculateMassOrMoleFluxes()
Calculate the mass or mole fluxes over all surfaces.
Definition: fluxoversurface.hh:225
FluxOverSurface(const GridVariables &gridVariables, const Sol &sol)
The constructor.
Definition: fluxoversurface.hh:148
static SurfaceGeometryType makeSurface(const std::vector< Dune::FieldVector< Scalar, 3 > > &corners)
Creates a geometrical surface object for (3D).
Definition: fluxoversurface.hh:217
void calculateFluxes(const FluxType &fluxType)
Calculate the fluxes over all surfaces for a given flux type.
Definition: fluxoversurface.hh:281
Detect if a point intersects a geometry.
Create Dune geometries from user-specified points.