From 4056c771dc7b9876653e77938fb2933f943811b4 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 12 Jul 2016 10:53:48 +0200
Subject: [PATCH] [fvGeometry] introduce global-local-bind concept
The global object can create an empty temporary local object that
can be bound to an element/stencil. This is closer to the
current dumux-stable implementation and still allows to switch
between global caching (for smaller problems where memory is not an issue)
and local caching (variables are precomputed for each element assemble temporarily)
Works for cc and box (in 2d). Box needs some performance improvement,
probably when creating the box geometries. But also other parts
might be possbile to optimize.
---
dumux/common/pointsource.hh | 8 +-
dumux/discretization/box/boxgeometryhelper.hh | 444 +++++++++
dumux/discretization/box/darcyslaw.hh | 77 +-
.../box/fluxvariablescachevector.hh | 126 ++-
dumux/discretization/box/fvelementgeometry.hh | 373 +++++++
.../box/fvelementgeometryvector.hh | 919 ------------------
dumux/discretization/box/globalfvgeometry.hh | 383 ++++++++
dumux/discretization/box/stencils.hh | 8 +-
dumux/discretization/box/subcontrolvolume.hh | 21 +-
.../box/subcontrolvolumeface.hh | 14 +-
.../box/volumevariablesvector.hh | 344 ++-----
.../cellcentered/fluxvariablescachevector.hh | 57 +-
dumux/discretization/cellcentered/stencils.hh | 26 +-
.../cellcentered/tpfa/CMakeLists.txt | 3 +-
.../cellcentered/tpfa/darcyslaw.hh | 58 +-
.../cellcentered/tpfa/fickslaw.hh | 22 +-
.../cellcentered/tpfa/fvelementgeometry.hh | 360 +++++++
...tgeometryvector.hh => globalfvgeometry.hh} | 248 ++---
.../cellcentered/tpfa/subcontrolvolume.hh | 11 +-
.../cellcentered/tpfa/subcontrolvolumeface.hh | 18 +-
.../cellcentered/volumevariablesvector.hh | 44 +-
dumux/discretization/fluxvariablesbase.hh | 11 +-
dumux/discretization/fvelementgeometry.hh | 207 ----
dumux/discretization/scvandscvfiterators.hh | 112 +++
dumux/discretization/subcontrolvolumebase.hh | 141 +--
.../subcontrolvolumefacebase.hh | 23 +-
dumux/geomechanics/elastic/model.hh | 1 +
dumux/implicit/assembler.hh | 3 +
dumux/implicit/box/elementboundarytypes.hh | 4 +-
dumux/implicit/box/localjacobian.hh | 81 +-
dumux/implicit/box/localresidual.hh | 126 +--
dumux/implicit/box/propertydefaults.hh | 13 +-
.../cellcentered/elementboundarytypes.hh | 4 +-
dumux/implicit/cellcentered/localjacobian.hh | 63 +-
dumux/implicit/cellcentered/localresidual.hh | 198 ++--
.../cellcentered/tpfa/propertydefaults.hh | 10 +-
dumux/implicit/localresidual.hh | 316 +++---
dumux/implicit/model.hh | 35 +-
dumux/implicit/problem.hh | 4 +-
dumux/implicit/properties.hh | 6 +-
dumux/implicit/propertydefaults.hh | 6 +-
dumux/porousmediumflow/1p/implicit/model.hh | 13 +-
dumux/porousmediumflow/2p/implicit/model.hh | 9 +-
dumux/porousmediumflow/3p/implicit/model.hh | 9 +-
.../3p3c/implicit/localresidual.hh | 11 +-
dumux/porousmediumflow/3p3c/implicit/model.hh | 18 +-
.../compositional/primaryvariableswitch.hh | 15 +-
.../immiscible/localresidual.hh | 31 +-
.../implicit/fluxvariables.hh | 44 +-
.../implicit/fluxvariablescache.hh | 29 +-
.../elastic/elasticmatrixproblem.hh | 2 +-
test/implicit/test_fvelementgeometry.cc | 6 +-
.../3p3c/implicit/infiltrationproblem.hh | 6 +-
53 files changed, 2697 insertions(+), 2424 deletions(-)
create mode 100644 dumux/discretization/box/boxgeometryhelper.hh
create mode 100644 dumux/discretization/box/fvelementgeometry.hh
delete mode 100644 dumux/discretization/box/fvelementgeometryvector.hh
create mode 100644 dumux/discretization/box/globalfvgeometry.hh
create mode 100644 dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
rename dumux/discretization/cellcentered/tpfa/{fvelementgeometryvector.hh => globalfvgeometry.hh} (59%)
delete mode 100644 dumux/discretization/fvelementgeometry.hh
create mode 100644 dumux/discretization/scvandscvfiterators.hh
diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 426567da3a..f636d72f91 100644
--- a/dumux/common/pointsource.hh
+++ b/dumux/common/pointsource.hh
@@ -314,14 +314,16 @@ public:
// check in which subcontrolvolume(s) we are
// TODO mapper/problem in bboxtree would allow to make this much better
const auto element = boundingBoxTree.entity(eIdx);
- auto fvGeometry = problem.model().fvGeometries(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
+
const auto globalPos = source.position();
// loop over all sub control volumes and check if the point source is inside
std::vector<unsigned int> scvIndices;
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
if (BoundingBoxTreeHelper<dimworld>::pointInGeometry(scv.geometry(), globalPos))
- scvIndices.push_back(scv.indexInElement());
+ scvIndices.push_back(scv.index());
}
// for all scvs that where tested positiv add the point sources
// to the element/scv to point source map
diff --git a/dumux/discretization/box/boxgeometryhelper.hh b/dumux/discretization/box/boxgeometryhelper.hh
new file mode 100644
index 0000000000..75e879c259
--- /dev/null
+++ b/dumux/discretization/box/boxgeometryhelper.hh
@@ -0,0 +1,444 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Helper class constructing the dual grid finite volume geometries
+ * for the box discretizazion method
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
+#define DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+//! Create sub control volumes and sub control volume face geometries
+template<class GridView, int dim>
+class BoxGeometryHelper;
+
+//! A class to create sub control volume and sub control volume face geometries per element
+template <class GridView>
+class BoxGeometryHelper<GridView, 1>
+{
+private:
+ using Scalar = typename GridView::ctype;
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
+ using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
+
+ using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+ using CornerList = std::vector<GlobalPosition>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
+
+ using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
+
+public:
+ using ScvGeometryVector = std::vector<ScvGeometry>;
+ using ScvfGeometryVector = std::vector<ScvfGeometry>;
+
+ //! get sub control volume geometries from element
+ void getScvAndScvfGeometries(const typename Element::Geometry& geometry,
+ ScvGeometryVector& scvGeometries,
+ ScvfGeometryVector& scvfGeometries)
+ {
+ scvGeometries.reserve(2);
+ scvfGeometries.reserve(1);
+
+ // the sub control volumes
+ scvGeometries.emplace_back(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.corner(0), geometry.center()}));
+ scvGeometries.emplace_back(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.center(), geometry.corner(1)}));
+
+ // the sub control volume faces
+ scvfGeometries.emplace_back(Dune::GeometryType(0), std::vector<GlobalPosition>({geometry.center()}));
+ }
+
+ //! get sub control volume geometries from element of dimension 1
+ ScvfGeometryVector getBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
+ {
+ return {ScvfGeometry(Dune::GeometryType(0), std::vector<GlobalPosition>({geometry.center()}))};
+ }
+
+ //! get scvf normal vector
+ static GlobalPosition normal(const typename Element::Geometry& geometry,
+ const ScvfGeometry& scvfGeometry)
+ {
+ GlobalPosition normal = geometry.corner(1) - geometry.corner(0);
+ normal /= normal.two_norm();
+ return normal;
+ }
+};
+
+//! A class to create sub control volume and sub control volume face geometries per element
+template <class GridView>
+class BoxGeometryHelper<GridView, 2>
+{
+private:
+ using Scalar = typename GridView::ctype;
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
+ using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
+
+ using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
+
+ using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
+
+public:
+
+ BoxGeometryHelper(const typename Element::Geometry& geometry)
+ {
+ // extract the corners of the sub control volumes
+ const auto& referenceElement = ReferenceElements::general(geometry.type());
+
+ // vertices
+ for (int i = 0; i < geometry.corners(); ++i)
+ v.emplace_back(geometry.corner(i));
+
+ // face midpoints
+ for (int i = 0; i < referenceElement.size(1); ++i)
+ f.emplace_back(geometry.global(referenceElement.position(i, 1)));
+
+ c = geometry.center();
+
+ corners = geometry.corners();
+ }
+
+ //! Create the sub control volume geometries
+ std::vector<ScvGeometry> createScvGeometries()
+ {
+ std::vector<ScvGeometry> scvGeometries;
+
+ // sub control volume geometries in 2D are always quadrilaterals
+ Dune::GeometryType scvGeometryType;
+ scvGeometryType.makeQuadrilateral();
+
+ // proceed according to number of corners
+ switch (corners)
+ {
+ case 3: // triangle
+ {
+ scvGeometries.reserve(3);
+
+ // the sub control volumes
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], f[0], f[1], c}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[0], v[1], c, f[2]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[1], c, v[2], f[2]}));
+
+ break;
+ }
+ case 4: // quadrilateral
+ {
+ scvGeometries.reserve(4);
+
+ // the sub control volumes
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], f[2], f[0], c}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[2], v[1], c, f[1]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[0], c, v[2], f[3]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({c, f[1], f[3], v[3]}));
+
+ break;
+ }
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
+ << " dimWorld=" << dimWorld
+ << " corners=" << corners);
+ }
+
+ return scvGeometries;
+ }
+
+
+ //! Create the sub control volume face geometries
+ std::vector<ScvfGeometry> createScvfGeometries()
+ {
+ std::vector<ScvfGeometry> scvfGeometries;
+
+ // sub control volume face geometries in 2D are always lines
+ Dune::GeometryType scvfGeometryType;
+ scvfGeometryType.makeLine();
+
+ // proceed according to number of corners
+ switch (corners)
+ {
+ case 3: // triangle
+ {
+ scvfGeometries.reserve(3);
+
+ // the sub control volume faces
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
+
+ break;
+ }
+ case 4: // quadrilateral
+ {
+ scvfGeometries.reserve(4);
+
+ // the sub control volume faces
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[3], c}));
+
+ break;
+ }
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
+ << " dimWorld=" << dimWorld
+ << " corners=" << corners);
+ }
+
+ return scvfGeometries;
+ }
+
+ //! Create the sub control volume face geometries on the boundary
+ static std::vector<ScvfGeometry> createBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
+ {
+ return {ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.corner(0), geometry.center()})),
+ ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.center(), geometry.corner(1)}))};
+ }
+
+ //! get scvf normal vector for dim == 2, dimworld == 3
+ template <int w = dimWorld>
+ static typename std::enable_if<w == 3, GlobalPosition>::type
+ normal(const typename Element::Geometry& geometry,
+ const ScvfGeometry& scvfGeometry)
+ {
+ const auto v1 = geometry.corner(1) - geometry.corner(0);
+ const auto v2 = geometry.corner(2) - geometry.corner(0);
+ const auto v3 = Dumux::crossProduct(v1, v2);
+ const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
+ GlobalPosition normal = Dumux::crossProduct(v3, t);
+ normal /= normal.two_norm();
+ return normal;
+ }
+
+ //! get scvf normal vector for dim == 2, dimworld == 2
+ template <int w = dimWorld>
+ static typename std::enable_if<w == 2, GlobalPosition>::type
+ normal(const typename Element::Geometry& geometry,
+ const ScvfGeometry& scvfGeometry)
+ {
+ const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
+ GlobalPosition normal({-t[1], t[0]});
+ normal /= normal.two_norm();
+ return normal;
+ }
+private:
+ std::vector<GlobalPosition> v; // vertices
+ std::vector<GlobalPosition> f; // face midpoints
+ GlobalPosition c; // element center
+ std::size_t corners; // number of element corners
+};
+
+//! A class to create sub control volume and sub control volume face geometries per element
+template <class GridView>
+class BoxGeometryHelper<GridView, 3>
+{
+private:
+ using Scalar = typename GridView::ctype;
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
+ using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
+
+ using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+ using CornerList = std::vector<GlobalPosition>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Intersection = typename GridView::Intersection;
+
+ using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
+ using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
+
+public:
+ using ScvGeometryVector = std::vector<ScvGeometry>;
+ using ScvfGeometryVector = std::vector<ScvfGeometry>;
+
+ //! get sub control volume geometries from element
+ void getScvAndScvfGeometries(const typename Element::Geometry& geometry,
+ ScvGeometryVector& scvGeometries,
+ ScvfGeometryVector& scvfGeometries)
+ {
+ // sub control volume geometries in 3D are always hexahedrons
+ Dune::GeometryType scvGeometryType; scvGeometryType.makeHexahedron();
+ // sub control volume face geometries in 3D are always quadrilaterals
+ Dune::GeometryType scvfGeometryType; scvfGeometryType.makeQuadrilateral();
+
+ // extract the corners of the sub control volumes
+ const auto& referenceElement = ReferenceElements::general(geometry.type());
+
+ // vertices
+ CornerList v;
+ for (int i = 0; i < geometry.corners(); ++i)
+ v.emplace_back(geometry.corner(i));
+
+ // edge midpoints
+ CornerList e;
+ for (int i = 0; i < referenceElement.size(dim-1); ++i)
+ e.emplace_back(geometry.global(referenceElement.position(i, dim-1)));
+
+ // face midpoints
+ CornerList f;
+ for (int i = 0; i < referenceElement.size(1); ++i)
+ f.emplace_back(geometry.global(referenceElement.position(i, 1)));
+
+ auto c = geometry.center();
+
+ // procees according to number of corners
+ // \todo prisms (corners == 6) and pyramids (corners == 5)
+ switch (geometry.corners())
+ {
+ case 4: // tetrahedron
+ {
+ scvGeometries.reserve(4);
+ scvfGeometries.reserve(6);
+
+ // sub control volumes
+ ScvGeometryVector scvGeometries(4);
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], e[0], e[1], f[0], e[3], f[1], f[2], c}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[1], e[2], e[0], f[0], f[3], e[4], c, f[1]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[2], e[1], e[2], f[0], e[5], f[2], f[3], c}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[3], e[3], e[5], f[2], e[4], f[1], f[3], c}));
+
+ // sub control volume faces
+ ScvfGeometryVector scvfGeometries(6);
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[0], f[0], f[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], e[1], c, f[2]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[2], f[0], f[3], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[2], e[3], c, f[1]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[3], c, e[4], f[1]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[5], f[2], f[3], c}));
+
+ break;
+ }
+ case 8: // hexahedron
+ {
+ scvGeometries.reserve(8);
+ scvfGeometries.reserve(12);
+
+ // sub control volumes
+ ScvGeometryVector scvGeometries(8);
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({v[0], e[6], e[4], f[4], e[0], f[2], f[0], c}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({e[6], v[1], f[4], e[5], f[2], e[1], c, f[1]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({e[4], f[4], v[2], e[7], f[0], c, e[2], f[3]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[4], e[5], e[7], v[3], c, f[1], f[3], e[3]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({e[0], f[2], f[0], c, v[4], e[10], e[8], f[5]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[2], e[1], c, f[1], e[10], v[5], f[5], e[9]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({f[0], c, e[2], f[3], e[8], f[5], v[6], e[11]}));
+ scvGeometries.emplace_back(scvGeometryType, std::vector<GlobalPosition>({c, f[1], f[3], e[3], f[5], e[9], e[11], v[7]}));
+
+ // sub control volume faces
+ ScvfGeometryVector scvfGeometries(12);
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[0], e[0], c, f[2]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[1], c, e[1], f[2]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[3], e[2], c, f[0]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[3], f[3], f[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[4], e[4], c, f[0]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[5], f[4], f[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[6], f[4], f[2], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({f[4], e[7], c, f[3]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({c, f[0], f[5], e[8]}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[9], f[1], f[5], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[10], f[2], f[5], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({e[11], f[5], f[3], c}));
+
+ break;
+ }
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Box scv geometries for " << geometry.type());
+ }
+ }
+
+ //! get sub control volume geometries from element
+ ScvfGeometryVector getBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
+ {
+ ScvfGeometryVector scvfGeometries;
+
+ // sub control volume face geometries in 3D are always quadrilaterals
+ Dune::GeometryType scvfGeometryType; scvfGeometryType.makeQuadrilateral();
+
+ // extract the corners of the sub control volumes
+ const auto& referenceElement = FaceReferenceElements::general(geometry.type());
+
+ // vertices
+ CornerList v;
+ for (int i = 0; i < geometry.corners(); ++i)
+ v.emplace_back(geometry.corner(i));
+
+ // edge midpoints
+ CornerList e;
+ for (int i = 0; i < referenceElement.size(1); ++i)
+ e.emplace_back(geometry.global(referenceElement.position(i, 1)));
+
+ // face midpoint
+ auto c = geometry.center();
+
+ // procees according to number of corners
+ switch (geometry.corners())
+ {
+ case 3: // triangle
+ {
+ scvfGeometries.reserve(3);
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[0], e[0], e[1], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[1], e[2], e[0], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[2], e[1], e[2], c}));
+
+ return scvfGeometries;
+ }
+ case 4: // quadrilateral
+ {
+ scvfGeometries.reserve(4);
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[0], e[2], e[0], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[1], e[1], e[2], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[2], e[0], e[3], c}));
+ scvfGeometries.emplace_back(scvfGeometryType, std::vector<GlobalPosition>({v[3], e[3], e[1], c}));
+
+ return scvfGeometries;
+ }
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Box scvf boundary geometries for " << geometry.type());
+ }
+ }
+
+ //! get scvf normal vector
+ static GlobalPosition normal(const typename Element::Geometry& geometry,
+ const ScvfGeometry& scvfGeometry)
+ {
+ const auto t1 = scvfGeometry.corner(1) - scvfGeometry.corner(0);
+ const auto t2 = scvfGeometry.corner(2) - scvfGeometry.corner(0);
+ GlobalPosition normal = Dumux::crossProduct(t1, t2);
+ normal /= normal.two_norm();
+ return normal;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/box/darcyslaw.hh b/dumux/discretization/box/darcyslaw.hh
index 61152459b0..093ba8ddf2 100644
--- a/dumux/discretization/box/darcyslaw.hh
+++ b/dumux/discretization/box/darcyslaw.hh
@@ -52,50 +52,47 @@ NEW_PROP_TAG(ProblemEnableGravity);
template <class TypeTag>
class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::Box>::type >
{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::IndexSet::IndexType IndexType;
- typedef typename GridView::ctype CoordScalar;
- typedef std::vector<IndexType> Stencil;
-
- enum { dim = GridView::dimension} ;
- enum { dimWorld = GridView::dimensionworld} ;
-
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
-
- typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> FeCache;
- typedef typename FeCache::FiniteElementType::Traits::LocalBasisType FeLocalBasis;
- typedef typename FluxVariablesCache::FaceData FaceData;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using CoordScalar = typename GridView::ctype;
+ using Stencil = std::vector<IndexType>;
+
+ enum { dim = GridView::dimension};
+ enum { dimWorld = GridView::dimensionworld};
+
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+ using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+ using FaceData = typename FluxVariablesCache::FaceData;
public:
static Scalar flux(const Problem& problem,
const Element& element,
- const SubControlVolumeFace& scvFace,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
const IndexType phaseIdx)
{
// get the precalculated local jacobian and shape values at the integration point
- const auto& faceData = problem.model().fluxVarsCache()[scvFace].faceData();
+ const auto& faceData = problem.model().fluxVarsCache()[scvf].faceData();
- const auto& fvGeometry = problem.model().fvGeometries(element);
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto extrusionFactor = problem.model().curVolVars(insideScv).extrusionFactor();
const auto K = problem.spatialParams().intrinsicPermeability(insideScv);
// evaluate gradP - rho*g at integration point
DimVector gradP(0.0);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
// the global shape function gradient
DimVector gradI;
- faceData.jacInvT.mv(faceData.localJacobian[scv.indexInElement()][0], gradI);
+ faceData.jacInvT.mv(faceData.shapeJacobian[scv.index()][0], gradI);
gradI *= problem.model().curVolVars(scv).pressure(phaseIdx);
gradP += gradI;
@@ -103,11 +100,11 @@ public:
if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
{
// gravitational acceleration
- DimVector g(problem.gravityAtPos(scvFace.center()));
+ DimVector g(problem.gravityAtPos(scvf.center()));
// interpolate the density at the IP
const auto& insideVolVars = problem.model().curVolVars(insideScv);
- const auto& outsideScv = problem.model().fvGeometries().subControlVolume(scvFace.outsideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& outsideVolVars = problem.model().curVolVars(outsideScv);
Scalar rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
@@ -120,7 +117,7 @@ public:
// apply the permeability and return the flux
auto KGradP = applyPermeability(K, gradP);
- return -1.0*(KGradP*scvFace.unitOuterNormal())*scvFace.area()*extrusionFactor;
+ return -1.0*(KGradP*scvf.unitOuterNormal())*scvf.area()*extrusionFactor;
}
static DimVector applyPermeability(const DimWorldMatrix& K, const DimVector& gradI)
@@ -139,22 +136,30 @@ public:
}
// This is for compatibility with the cc methods. The flux stencil info is obsolete for the box method.
- static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+ static Stencil stencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
return Stencil(0);
}
- static FaceData calculateFaceData(const Problem& problem, const Element& element, const typename Element::Geometry& geometry, const FeLocalBasis& localBasis, const SubControlVolumeFace& scvFace)
+ static FaceData calculateFaceData(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
FaceData faceData;
+ const auto geometry = element.geometry();
+ const auto& localBasis = fvGeometry.feLocalBasis();
// evaluate shape functions and gradients at the integration point
- const auto ipLocal = geometry.local(scvFace.center());
+ const auto ipLocal = geometry.local(scvf.center());
faceData.jacInvT = geometry.jacobianInverseTransposed(ipLocal);
- localBasis.evaluateJacobian(ipLocal, faceData.localJacobian);
- localBasis.evaluateFunction(ipLocal, faceData.shapeValues);
+ localBasis.evaluateJacobian(ipLocal, faceData.shapeJacobian);
+ localBasis.evaluateFunction(ipLocal, faceData.shapeValue);
- return std::move(faceData);
+ return faceData;
}
};
diff --git a/dumux/discretization/box/fluxvariablescachevector.hh b/dumux/discretization/box/fluxvariablescachevector.hh
index b62bc5e92c..65168271c8 100644
--- a/dumux/discretization/box/fluxvariablescachevector.hh
+++ b/dumux/discretization/box/fluxvariablescachevector.hh
@@ -28,12 +28,22 @@
namespace Dumux
{
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for box models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFluxVariablesCache>
+class BoxFluxVariablesCacheVector
+{};
+
/*!
* \ingroup ImplicitModel
* \brief Base class for the flux variables cache vector, we store one cache per face
*/
template<class TypeTag>
-class BoxFluxVariablesCacheVector
+class BoxFluxVariablesCacheVector<TypeTag, true>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -41,87 +51,107 @@ class BoxFluxVariablesCacheVector
using Element = typename GridView::template Codim<0>::Entity;
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
public:
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ void update(Problem& problem)
{
problemPtr_ = &problem;
- fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
+ fluxVarsCache_.resize(problem.gridView.size(0));
for (const auto& element : elements(problem.gridView()))
{
+ eIdx_ = problem.elementMapper().index(element);
// bind the geometries and volume variables to the element (all the elements in stencil)
- problem.model().fvGeometries_().bind(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometries());
+ fvGeometry.bind(element);
- const auto& elementGeometry = element.geometry();
- const auto& localBasis = problem.model().fvGeometries().feLocalBasis(elementGeometry.type());
- const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scvf : scvfs(fvGeometry))
+ const auto& localBasis = fvGeometry.feLocalBasis();
+ fluxVarsCache_[eIdx_].resize(fvGeometry.numScvf());
+ for (auto&& scvf : scvfs(fvGeometry))
{
- (*this)[scvf].update(problem, element, elementGeometry, localBasis, scvf);
+ (*this)[scvf].update(problem, element, localBasis, scvf);
}
}
}
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ // Function is called by the BoxLocalJacobian prior to flux calculations on the element.
+ // We assume the FVGeometries to be bound at this point
+ void bind(const Element& element, const FVElementGeometry& fvGeometry)
+ {
+ bindElement(element, fvGeometry);
+ }
+
+ void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
+ { eIdx_ = problem_().elementMapper().index(element); }
+
+ // access operator
+ const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
+ { return fluxVarsCache_[eIdx_][scvf.index()]; }
+
+ // access operator
+ FluxVariablesCache& operator [](const SubControlVolumeFace& scvf)
+ { return fluxVarsCache_[eIdx_][scvf.index()]; }
+
+private:
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ // currently bound element
+ IndexType eIdx_;
+ const Problem* problemPtr_;
+ std::vector<std::vector<FluxVariablesCache>> fluxVarsCache_;
+};
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the flux variables cache vector, we store one cache per face
+ */
+template<class TypeTag>
+class BoxFluxVariablesCacheVector<TypeTag, false>
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+public:
+
+ void update(Problem& problem)
{
problemPtr_ = &problem;
}
// Function is called by the BoxLocalJacobian prior to flux calculations on the element.
// We assume the FVGeometries to be bound at this point
- void bind(const Element& element)
+ void bind(const Element& element, const FVElementGeometry& fvGeometry)
{
- bindElement(element);
+ bindElement(element, fvGeometry);
}
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
- bindElement(const Element& element)
+ void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
{
- const auto& fvGeometry = problem_().model().fvGeometries(element);
- const auto& elementGeometry = element.geometry();
- const auto& localBasis = problem_().model().fvGeometries().feLocalBasis(elementGeometry.type());
-
// temporary resizing of the cache
- const auto numScvf = fvGeometry.numScvf();
- fluxVarsCache_.resize(numScvf);
- IndexType localScvfIdx = 0;
- for (const auto& scvf : scvfs(fvGeometry))
- fluxVarsCache_[localScvfIdx++].update(problem_(), element, elementGeometry, localBasis, scvf);
+ fluxVarsCache_.resize(fvGeometry.numScvf());
+ for (auto&& scvf : scvfs(fvGeometry))
+ (*this)[scvf].update(problem_(), element, fvGeometry, scvf);
}
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
- bindElement(const Element& element) {}
-
- // access operators in the case of caching
- template <typename T = TypeTag>
- const typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
- operator [](const SubControlVolumeFace& scvf) const
+ // access operator
+ const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
{ return fluxVarsCache_[scvf.index()]; }
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
- operator [](const SubControlVolumeFace& scvf)
+ // access operator
+ FluxVariablesCache& operator [](const SubControlVolumeFace& scvf)
{ return fluxVarsCache_[scvf.index()]; }
- // access operators in the case of no caching
- template <typename T = TypeTag>
- const typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
- operator [](const SubControlVolumeFace& scvf) const
- { return fluxVarsCache_[scvf.indexInElement()]; }
-
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache), FluxVariablesCache>::type&
- operator [](const SubControlVolumeFace& scvf)
- { return fluxVarsCache_[scvf.indexInElement()]; }
-
private:
const Problem& problem_() const
{ return *problemPtr_; }
+ // currently bound element
const Problem* problemPtr_;
std::vector<FluxVariablesCache> fluxVarsCache_;
};
diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
new file mode 100644
index 0000000000..5d843b4809
--- /dev/null
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -0,0 +1,373 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the local finite volume geometry for box models
+ * This builds up the sub control volumes and sub control volume faces
+ * for an element.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_FV_ELEMENT_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_BOX_FV_ELEMENT_GEOMETRY_HH
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dumux/discretization/scvandscvfiterators.hh>
+#include <dumux/discretization/box/boxgeometryhelper.hh>
+#include <dumux/implicit/box/properties.hh>
+
+namespace Dumux
+{
+
+//! forward declaration of the global finite volume geometry
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class BoxGlobalFVGeometry;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for box models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class BoxFVElementGeometry
+{};
+
+//! specialization in case the FVElementGeometries are stored
+template<class TypeTag>
+class BoxFVElementGeometry<TypeTag, true>
+{
+ using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using CoordScalar = typename GridView::ctype;
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+
+ using ScvIterator = ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
+ using ScvfIterator = ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
+
+public:
+ //! Constructor
+ BoxFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
+ : globalFvGeometryPtr_(&globalFvGeometry) {}
+
+ //! Get a sub control volume with a local scv index
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ return globalFvGeometry().scv(scvIdx, eIdx_);
+ }
+
+ //! Get a sub control volume face with a local scvf index
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ return globalFvGeometry().scvf(scvfIdx, eIdx_);
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
+ scvs(const BoxFVElementGeometry& fvGeometry)
+ {
+ const auto& g = fvGeometry.globalFvGeometry();
+ using Iter = typename std::vector<SubControlVolume>::const_iterator;
+ return Dune::IteratorRange<Iter>(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
+ scvfs(const BoxFVElementGeometry& fvGeometry)
+ {
+ const auto& g = fvGeometry.globalFvGeometry();
+ using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<Iter>(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end());
+ }
+
+ //! Get a local finite element basis
+ const FeLocalBasis& feLocalBasis() const
+ {
+ return globalFvGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+ }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return globalFvGeometry().scvs(eIdx_).size();
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return globalFvGeometry().scvfs(eIdx_).size();
+ }
+
+ //! this function is for compatibility reasons with cc methods
+ //! The box stencil is always element-local so bind and bindElement
+ //! are identical.
+ void bind(const Element& element)
+ {
+ this->bindElement(element);
+ }
+
+ //! Binding of an element, has to be called before using the fvgeometries
+ //! Prepares all the volume variables within the element
+ //! For compatibility reasons with the FVGeometry cache being disabled
+ void bindElement(const Element& element)
+ {
+ elementPtr_ = &element;
+ eIdx_ = globalFvGeometry().problem_().elementMapper().index(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ const GlobalFVGeometry& globalFvGeometry() const
+ { return *globalFvGeometryPtr_; }
+
+private:
+ const Element* elementPtr_;
+ const GlobalFVGeometry* globalFvGeometryPtr_;
+
+ IndexType eIdx_;
+};
+
+//! specialization in case the FVElementGeometries are not stored
+template<class TypeTag>
+class BoxFVElementGeometry<TypeTag, false>
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using CoordScalar = typename GridView::ctype;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+
+ using GeometryHelper = BoxGeometryHelper<GridView, dim>;
+
+public:
+ //! Constructor
+ BoxFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
+ : globalFvGeometryPtr_(&globalFvGeometry) {}
+
+ //! Get a sub control volume with a local scv index
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ return scvs_[scvIdx];
+ }
+
+ //! Get a sub control volume face with a local scvf index
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ return scvfs_[scvfIdx];
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
+ scvs(const BoxFVElementGeometry& fvGeometry)
+ {
+ using Iter = typename std::vector<SubControlVolume>::const_iterator;
+ return Dune::IteratorRange<Iter>(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element.
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
+ scvfs(const BoxFVElementGeometry& fvGeometry)
+ {
+ using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<Iter>(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end());
+ }
+
+ //! Get a local finite element basis
+ const FeLocalBasis& feLocalBasis() const
+ {
+ return globalFvGeometry().feCache().get(elementPtr_->geometry().type()).localBasis();
+ }
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return scvs_.size();
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return scvfs_.size();
+ }
+
+ //! this function is for compatibility reasons with cc methods
+ //! The box stencil is always element-local so bind and bindElement
+ //! are identical.
+ void bind(const Element& element)
+ {
+ this->bindElement(element);
+ }
+
+ //! Binding of an element, has to be called before using the fvgeometries
+ //! Prepares all the volume variables within the element
+ //! For compatibility reasons with the FVGeometry cache being disabled
+ void bindElement(const Element& element)
+ {
+ elementPtr_ = &element;
+ eIdx_ = globalFvGeometry().problem_().elementMapper().index(element);
+ makeElementGeometries(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ const GlobalFVGeometry& globalFvGeometry() const
+ { return *globalFvGeometryPtr_; }
+
+private:
+
+ void makeElementGeometries(const Element& element)
+ {
+ auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
+
+ // get the element geometry
+ auto elementGeometry = element.geometry();
+ const auto& referenceElement = ReferenceElements::general(elementGeometry.type());
+
+ // get the sub control volume geometries of this element
+ GeometryHelper geometryHelper(elementGeometry);
+ auto scvGeometries = geometryHelper.createScvGeometries();
+ auto scvfGeometries = geometryHelper.createScvfGeometries();
+
+ // construct the sub control volumes
+ IndexType scvLocalIdx = 0;
+ for (auto&& scvGeometry : scvGeometries)
+ {
+ // get asssociated dof index
+ auto dofIdxGlobal = globalFvGeometry().problem_().vertexMapper().subIndex(element, scvLocalIdx, dim);
+
+ // add scv to the local container
+ scvs_.emplace_back(std::move(scvGeometry),
+ scvLocalIdx,
+ eIdx,
+ dofIdxGlobal);
+ // increment local counter
+ scvLocalIdx++;
+ }
+
+ // construct the sub control volume faces
+ IndexType scvfLocalIdx = 0;
+ for (auto&& scvfGeometry : scvfGeometries)
+ {
+ // find the local scv indices this scvf is connsected to
+ std::vector<IndexType> localScvIndices({static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+ static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+
+ // compute the scvf normal unit outer normal
+ auto normal = geometryHelper.normal(elementGeometry, scvfGeometry);
+ const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
+ const auto s = v*normal;
+ if (std::signbit(s))
+ normal *= -1;
+
+ scvfs_.emplace_back(std::move(scvfGeometry),
+ normal,
+ scvfLocalIdx,
+ localScvIndices,
+ false);
+
+ // increment local counter
+ scvfLocalIdx++;
+ }
+
+ // construct the sub control volume faces on the domain boundary
+ for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
+ {
+ if (intersection.boundary())
+ {
+ auto isGeometry = intersection.geometry();
+ auto boundaryScvfGeometries = geometryHelper.createBoundaryScvfGeometries(isGeometry);
+
+ IndexType boundaryFaceLocalIdx = 0;
+ for (auto&& scvfGeometry : boundaryScvfGeometries)
+ {
+ // find the scv this scvf is connected to
+ std::vector<IndexType> localScvIndices =
+ {static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1,
+ boundaryFaceLocalIdx++, dim))};
+
+ scvfs_.emplace_back(std::move(scvfGeometry),
+ intersection.centerUnitOuterNormal(),
+ scvfLocalIdx,
+ localScvIndices,
+ true);
+
+ // increment local counter
+ scvfLocalIdx++;
+ }
+ }
+ }
+ }
+
+ //! The bound element
+ const Element* elementPtr_;
+ IndexType eIdx_;
+
+ //! The global geometry this is a restriction of
+ const GlobalFVGeometry* globalFvGeometryPtr_;
+
+ //! vectors to store the geometries locally after binding an element
+ std::vector<SubControlVolume> scvs_;
+ std::vector<SubControlVolumeFace> scvfs_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/box/fvelementgeometryvector.hh b/dumux/discretization/box/fvelementgeometryvector.hh
deleted file mode 100644
index 48c11d229f..0000000000
--- a/dumux/discretization/box/fvelementgeometryvector.hh
+++ /dev/null
@@ -1,919 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Base class for the finite volume geometry vector for box models
- * This builds up the sub control volumes and sub control volume faces
- * for each element.
- */
-#ifndef DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
-#define DUMUX_DISCRETIZATION_BOX_FV_GEOMETRY_VECTOR_HH
-
-#include <dune/geometry/multilineargeometry.hh>
-#include <dune/geometry/referenceelements.hh>
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
-
-#include <dumux/implicit/box/properties.hh>
-#include <dumux/common/elementmap.hh>
-#include <dumux/common/math.hh>
-
-namespace Dumux
-{
-
-//! A class to create sub control volume and sub control volume face geometries per element
-template <class GridView>
-class BoxGeometryHelper
-{
-private:
- using Scalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
- using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
- using ScvfGeometry = Dune::MultiLinearGeometry<Scalar, dim-1, dimWorld>;
-
- using ScvGeometryVector = std::vector<std::shared_ptr<ScvGeometry>>;
- using ScvfGeometryVector = std::vector<std::shared_ptr<ScvfGeometry>>;
- using GeometryPair = std::pair<ScvGeometryVector, ScvfGeometryVector>;
-
- using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
- using CornerList = std::vector<GlobalPosition>;
- using Element = typename GridView::template Codim<0>::Entity;
- using Intersection = typename GridView::Intersection;
-
- using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
- using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
-
-public:
- //! get sub control volume geometries from element of dimension 1
- template <int d = dim>
- static typename std::enable_if<d == 1, GeometryPair>::type
- getScvAndScvfGeometries(const typename Element::Geometry& geometry)
- {
- // the sub control volumes
- ScvGeometryVector scvGeometries(2);
- scvGeometries[0] = std::make_shared<ScvGeometry>(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.corner(0), geometry.center()}));
- scvGeometries[1] = std::make_shared<ScvGeometry>(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.center(), geometry.corner(1)}));
-
- // the sub control volume faces
- ScvfGeometryVector scvfGeometries(1);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(Dune::GeometryType(0), std::vector<GlobalPosition>({geometry.center()}));
-
- return std::make_pair(scvGeometries, scvfGeometries);
- }
-
- //! get sub control volume geometries from element of dimension 2
- template <int d = dim>
- static typename std::enable_if<d == 2, GeometryPair>::type
- getScvAndScvfGeometries(const typename Element::Geometry& geometry)
- {
- // sub control volume geometries in 2D are always quadrilaterals
- Dune::GeometryType scvGeometryType; scvGeometryType.makeQuadrilateral();
- // sub control volume face geometries in 2D are always lines
- Dune::GeometryType scvfGeometryType; scvfGeometryType.makeLine();
-
- // extract the corners of the sub control volumes
- const auto& referenceElement = ReferenceElements::general(geometry.type());
-
- // vertices
- CornerList v;
- for (int i = 0; i < geometry.corners(); ++i)
- v.emplace_back(geometry.corner(i));
-
- // face midpoints
- CornerList f;
- for (int i = 0; i < referenceElement.size(1); ++i)
- f.emplace_back(geometry.global(referenceElement.position(i, 1)));
-
- auto c = geometry.center();
-
- // proceed according to number of corners
- switch (geometry.corners())
- {
- case 3: // triangle
- {
- // the sub control volumes
- ScvGeometryVector scvGeometries(3);
- scvGeometries[0] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[0], f[0], f[1], c}));
- scvGeometries[1] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[0], v[1], c, f[2]}));
- scvGeometries[2] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[1], c, v[2], f[2]}));
-
- // the sub control volume faces
- ScvfGeometryVector scvfGeometries(3);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
-
- return std::make_pair(scvGeometries, scvfGeometries);
- }
- case 4: // quadrilateral
- {
- // the sub control volumes
- ScvGeometryVector scvGeometries(4);
- scvGeometries[0] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[0], f[2], f[0], c}));
- scvGeometries[1] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[2], v[1], c, f[1]}));
- scvGeometries[2] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[0], c, v[2], f[3]}));
- scvGeometries[3] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({c, f[1], f[3], v[3]}));
-
- // the sub control volume faces
- ScvfGeometryVector scvfGeometries(4);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[0], c}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[1], c}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[2], c}));
- scvfGeometries[3] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[3], c}));
-
- return std::make_pair(scvGeometries, scvfGeometries);
- }
- default:
- DUNE_THROW(Dune::NotImplemented, "Box scv geometries for " << geometry.type());
- }
- }
-
- //! get sub control volume geometries from element of dimension 3
- template <int d = dim>
- static typename std::enable_if<d == 3, GeometryPair>::type
- getScvAndScvfGeometries(const typename Element::Geometry& geometry)
- {
- // sub control volume geometries in 3D are always hexahedrons
- Dune::GeometryType scvGeometryType; scvGeometryType.makeHexahedron();
- // sub control volume face geometries in 3D are always quadrilaterals
- Dune::GeometryType scvfGeometryType; scvfGeometryType.makeQuadrilateral();
-
- // extract the corners of the sub control volumes
- const auto& referenceElement = ReferenceElements::general(geometry.type());
-
- // vertices
- CornerList v;
- for (int i = 0; i < geometry.corners(); ++i)
- v.emplace_back(geometry.corner(i));
-
- // edge midpoints
- CornerList e;
- for (int i = 0; i < referenceElement.size(dim-1); ++i)
- e.emplace_back(geometry.global(referenceElement.position(i, dim-1)));
-
- // face midpoints
- CornerList f;
- for (int i = 0; i < referenceElement.size(1); ++i)
- f.emplace_back(geometry.global(referenceElement.position(i, 1)));
-
- auto c = geometry.center();
-
- // procees according to number of corners
- // \todo prisms (corners == 6) and pyramids (corners == 5)
- switch (geometry.corners())
- {
- case 4: // tetrahedron
- {
- // sub control volumes
- ScvGeometryVector scvGeometries(4);
- scvGeometries[0] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[0], e[0], e[1], f[0], e[3], f[1], f[2], c}));
- scvGeometries[1] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[1], e[2], e[0], f[0], f[3], e[4], c, f[1]}));
- scvGeometries[2] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[2], e[1], e[2], f[0], e[5], f[2], f[3], c}));
- scvGeometries[3] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[3], e[3], e[5], f[2], e[4], f[1], f[3], c}));
-
- // sub control volume faces
- ScvfGeometryVector scvfGeometries(6);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[0], f[0], f[1], c}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[0], e[1], c, f[2]}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[2], f[0], f[3], c}));
- scvfGeometries[3] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[2], e[3], c, f[1]}));
- scvfGeometries[4] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[3], c, e[4], f[1]}));
- scvfGeometries[5] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[5], f[2], f[3], c}));
-
- return std::make_pair(scvGeometries, scvfGeometries);
- }
- case 8: // hexahedron
- {
- // sub control volumes
- ScvGeometryVector scvGeometries(8);
- scvGeometries[0] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({v[0], e[6], e[4], f[4], e[0], f[2], f[0], c}));
- scvGeometries[1] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({e[6], v[1], f[4], e[5], f[2], e[1], c, f[1]}));
- scvGeometries[2] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({e[4], f[4], v[2], e[7], f[0], c, e[2], f[3]}));
- scvGeometries[3] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[4], e[5], e[7], v[3], c, f[1], f[3], e[3]}));
- scvGeometries[4] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({e[0], f[2], f[0], c, v[4], e[10], e[8], f[5]}));
- scvGeometries[5] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[2], e[1], c, f[1], e[10], v[5], f[5], e[9]}));
- scvGeometries[6] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({f[0], c, e[2], f[3], e[8], f[5], v[6], e[11]}));
- scvGeometries[7] = std::make_shared<ScvGeometry>(scvGeometryType, std::vector<GlobalPosition>({c, f[1], f[3], e[3], f[5], e[9], e[11], v[7]}));
-
- // sub control volume faces
- ScvfGeometryVector scvfGeometries(12);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[0], e[0], c, f[2]}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[1], c, e[1], f[2]}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[3], e[2], c, f[0]}));
- scvfGeometries[3] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[3], f[3], f[1], c}));
- scvfGeometries[4] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[4], e[4], c, f[0]}));
- scvfGeometries[5] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[5], f[4], f[1], c}));
- scvfGeometries[6] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[6], f[4], f[2], c}));
- scvfGeometries[7] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({f[4], e[7], c, f[3]}));
- scvfGeometries[8] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({c, f[0], f[5], e[8]}));
- scvfGeometries[9] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[9], f[1], f[5], c}));
- scvfGeometries[10] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[10], f[2], f[5], c}));
- scvfGeometries[11] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({e[11], f[5], f[3], c}));
-
- return std::make_pair(scvGeometries, scvfGeometries);
- }
- default:
- DUNE_THROW(Dune::NotImplemented, "Box scv geometries for " << geometry.type());
- }
- }
-
- //! get sub control volume geometries from element of dimension 1
- template <int d = dim>
- static typename std::enable_if<d == 1, ScvfGeometryVector>::type
- getBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
- {
- return {std::make_shared<ScvfGeometry>(ScvfGeometry(Dune::GeometryType(0), std::vector<GlobalPosition>({geometry.center()})))};
- }
-
- //! get sub control volume geometries from element of dimension 2
- template <int d = dim>
- static typename std::enable_if<d == 2, ScvfGeometryVector>::type
- getBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
- {
- return {std::make_shared<ScvfGeometry>(ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.corner(0), geometry.center()}))),
- std::make_shared<ScvfGeometry>(ScvfGeometry(Dune::GeometryType(1), std::vector<GlobalPosition>({geometry.center(), geometry.corner(1)})))};
- }
-
- //! get sub control volume geometries from element of dimension 3
- template <int d = dim>
- static typename std::enable_if<d == 3, ScvfGeometryVector>::type
- getBoundaryScvfGeometries(const typename Intersection::Geometry& geometry)
- {
- // sub control volume face geometries in 3D are always quadrilaterals
- Dune::GeometryType scvfGeometryType; scvfGeometryType.makeQuadrilateral();
-
- // extract the corners of the sub control volumes
- const auto& referenceElement = FaceReferenceElements::general(geometry.type());
-
- // vertices
- CornerList v;
- for (int i = 0; i < geometry.corners(); ++i)
- v.emplace_back(geometry.corner(i));
-
- // edge midpoints
- CornerList e;
- for (int i = 0; i < referenceElement.size(1); ++i)
- e.emplace_back(geometry.global(referenceElement.position(i, 1)));
-
- // face midpoint
- auto c = geometry.center();
-
- // procees according to number of corners
- switch (geometry.corners())
- {
- case 3: // triangle
- {
- ScvfGeometryVector scvfGeometries(3);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[0], e[0], e[1], c}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[1], e[2], e[0], c}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[2], e[1], e[2], c}));
-
- return scvfGeometries;
- }
- case 4: // quadrilateral
- {
- ScvfGeometryVector scvfGeometries(4);
- scvfGeometries[0] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[0], e[2], e[0], c}));
- scvfGeometries[1] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[1], e[1], e[2], c}));
- scvfGeometries[2] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[2], e[0], e[3], c}));
- scvfGeometries[3] = std::make_shared<ScvfGeometry>(scvfGeometryType, std::vector<GlobalPosition>({v[3], e[3], e[1], c}));
-
- return scvfGeometries;
- }
- default:
- DUNE_THROW(Dune::NotImplemented, "Box scvf boundary geometries for " << geometry.type());
- }
- }
-
- //! get scvf normal vector for dim == 3, dimworld == 3
- template <int d = dim, int w = dimWorld>
- static typename std::enable_if<d == 3 && w == 3, GlobalPosition>::type
- normal(const typename Element::Geometry& geometry, const ScvfGeometry& scvfGeometry)
- {
- const auto t1 = scvfGeometry.corner(1) - scvfGeometry.corner(0);
- const auto t2 = scvfGeometry.corner(2) - scvfGeometry.corner(0);
- GlobalPosition normal = Dumux::crossProduct(t1, t2);
- normal /= normal.two_norm();
- return normal;
- }
-
- //! get scvf normal vector for dim == 2, dimworld == 3
- template <int d = dim, int w = dimWorld>
- static typename std::enable_if<d == 2 && w == 3, GlobalPosition>::type
- normal(const typename Element::Geometry& geometry, const ScvfGeometry& scvfGeometry)
- {
- const auto v1 = geometry.corner(1) - geometry.corner(0);
- const auto v2 = geometry.corner(2) - geometry.corner(0);
- const auto v3 = Dumux::crossProduct(v1, v2);
- const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
- GlobalPosition normal = Dumux::crossProduct(v3, t);
- normal /= normal.two_norm();
- return normal;
- }
-
- //! get scvf normal vector for dim == 2, dimworld == 2
- template <int d = dim, int w = dimWorld>
- static typename std::enable_if<d == 2 && w == 2, GlobalPosition>::type
- normal(const typename Element::Geometry& geometry, const ScvfGeometry& scvfGeometry)
- {
- const auto t = scvfGeometry.corner(1) - scvfGeometry.corner(0);
- GlobalPosition normal({-t[1], t[0]});
- normal /= normal.two_norm();
- return normal;
- }
-
- //! get scvf normal vector for dim == 2, dimworld == 2
- template <int d = dim>
- static typename std::enable_if<d == 1, GlobalPosition>::type
- normal(const typename Element::Geometry& geometry, const ScvfGeometry& scvfGeometry)
- {
- GlobalPosition normal = geometry.corner(1) - geometry.corner(0);
- normal /= normal.two_norm();
- return normal;
- }
-};
-
-/*!
- * \ingroup ImplicitModel
- * \brief Base class for the finite volume geometry vector for box models
- * This builds up the sub control volumes and sub control volume faces
- * for each element.
- */
-template<class TypeTag, bool EnableFVElementGeometryCache>
-class BoxFVElementGeometryVector
-{};
-
-// specialization in case the FVElementGeometries are stored
-template<class TypeTag>
-class BoxFVElementGeometryVector<TypeTag, true>
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using Element = typename GridView::template Codim<0>::Entity;
-
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using CoordScalar = typename GridView::ctype;
-
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
- using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
- using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
-public:
- //! Constructor
- BoxFVElementGeometryVector(const GridView gridView)
- : gridView_(gridView), elementMap_(gridView) {}
-
- /* \brief Get the finite volume geometry of an element
- * \note The finite volume geometry offers iterators over the sub control volumes
- * and the sub control volume faces of an element.
- */
- const FVElementGeometry& fvGeometry(const Element& element) const
- {
- return fvGeometries_[problem_().elementMapper().index(element)];
- }
-
- //! Get a sub control volume with a global scv index
- const SubControlVolume& subControlVolume(IndexType scvIdx) const
- {
- return *scvs_[globalScvIndices_[eIdxBound_][scvIdx]];
- }
-
- //! Get a sub control volume face with a global scvf index
- const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
- {
- return *scvfs_[globalScvfIndices_[eIdxBound_][scvfIdx]];
- }
-
- //! Get a local finite element basis
- const FeLocalBasis& feLocalBasis(const Dune::GeometryType& type) const
- {
- return feCache_.get(type).localBasis();
- }
-
- //! The total number of sub control volumes
- std::size_t numScv() const
- {
- return scvs_.size();
- }
-
- //! The total number of sun control volume faces
- std::size_t numScvf() const
- {
- return scvfs_.size();
- }
-
- //! The total number of boundary sub control volume faces
- //! For compatibility reasons with cc methods
- std::size_t numBoundaryScvf() const
- {
- return 0;
- }
-
- // Get an element from a sub control volume contained in it
- Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
-
- // Get an element from a global element index
- Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
-
- //! update all fvElementGeometries (do this again after grid adaption)
- void update(const Problem& problem)
- {
- problemPtr_ = &problem;
-
- scvs_.clear();
- scvfs_.clear();
- fvGeometries_.clear();
- elementMap_.clear();
-
- IndexType numElements = gridView_.size(0);
- IndexType numScvs = 0;
- IndexType numScvfs = 0;
- for (const auto& element : elements(gridView_))
- {
- numScvs += element.subEntities(dim);
- numScvfs += element.subEntities(dim-1);
-
- // add number of boundary faces
- for (const auto& is : intersections(gridView_, element))
- {
- if (is.boundary())
- numScvfs += is.geometry().corners();
- }
- }
-
- scvs_.reserve(numScvs);
- scvfs_.reserve(numScvfs);
- globalScvIndices_.resize(numElements);
- globalScvfIndices_.resize(numElements);
- elementMap_.resize(numElements);
-
- // Build the SCV and SCV faces
- IndexType scvIdx = 0;
- IndexType scvfIdx = 0;
- for (const auto& element : elements(gridView_))
- {
- // the element-wise index sets for finite volume geometry
- std::vector<IndexType> scvLocalIndexSet;
- std::vector<IndexType> scvfLocalIndexSet;
-
- // fill the element map with seeds
- auto eIdx = problem.elementMapper().index(element);
- elementMap_[eIdx] = element.seed();
-
- // get the element geometry
- auto elementGeometry = element.geometry();
- const auto& referenceElement = ReferenceElements::general(elementGeometry.type());
-
- // get the sub control volume geometries of this element
- auto boxGeometries = std::move(BoxGeometryHelper<GridView>::getScvAndScvfGeometries(elementGeometry));
- // define some aliases
- auto& scvGeometries = boxGeometries.first;
- auto& scvfGeometries = boxGeometries.second;
-
- // construct the sub control volumes
- globalScvIndices_[eIdx].resize(scvGeometries.size());
- scvLocalIndexSet.reserve(scvGeometries.size());
- IndexType scvLocalIdx = 0;
- for (auto&& scvGeometry : scvGeometries)
- {
- scvLocalIndexSet.push_back(scvLocalIdx);
- auto dofIdxGlobal = problem.vertexMapper().subIndex(element, scvLocalIdx, dim);
- scvs_.emplace_back(std::make_shared<SubControlVolume>(*scvGeometry,
- scvIdx,
- eIdx,
- scvLocalIdx,
- dofIdxGlobal));
-
- globalScvIndices_[eIdx][scvLocalIdx] = scvIdx;
-
- // increment local counters
- scvIdx++; scvLocalIdx++;
- }
-
- // construct the sub control volume faces
- globalScvfIndices_[eIdx].resize(scvfGeometries.size());
- scvfLocalIndexSet.reserve(scvfGeometries.size());
- IndexType scvfLocalIdx = 0;
- for (auto&& scvfGeometry : scvfGeometries)
- {
- // add this scvf to the element's scvf list
- scvfLocalIndexSet.push_back(scvfLocalIdx);
-
- // find the global and local scv indices this scvf is belonging to
- std::vector<IndexType> localScvIndices(2);
- localScvIndices[0] = referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim);
- localScvIndices[1] = referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim);
-
- // compute the scvf normal unit outer normal
- auto normal = std::move(BoxGeometryHelper<GridView>::normal(elementGeometry, *scvfGeometry));
- const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
- const auto s = v*normal;
- if (std::signbit(s))
- normal *= -1;
-
- scvfs_.emplace_back(std::make_shared<SubControlVolumeFace>(*scvfGeometry,
- scvfGeometry->center(),
- normal,
- scvfIdx,
- scvfLocalIdx,
- localScvIndices,
- false));
-
- globalScvfIndices_[eIdx][scvfLocalIdx] = scvfIdx;
-
- // increment local counters
- scvfIdx++; scvfLocalIdx++;
- }
-
- // construct the sub control volume faces on the domain boundary
- for (const auto& intersection : intersections(gridView_, element))
- {
- if (intersection.boundary())
- {
- auto isGeometry = intersection.geometry();
- auto boundaryScvfGeometries = std::move(BoxGeometryHelper<GridView>::getBoundaryScvfGeometries(isGeometry));
-
- IndexType isScvfLocalIdx = 0;
- for (auto&& scvfGeometry : boundaryScvfGeometries)
- {
- // add this scvf to the element's scvf list
- scvfLocalIndexSet.push_back(scvfLocalIdx);
-
- // find the scvs this scvf is belonging to
- std::vector<IndexType> localScvIndices = {static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim))};
-
- // get the unit outer normal through the intersection
- auto normal = intersection.unitOuterNormal(isGeometry.local(scvfGeometry->center()));
-
- scvfs_.emplace_back(std::make_shared<SubControlVolumeFace>(*scvfGeometry,
- scvfGeometry->center(),
- normal,
- scvfIdx,
- scvfLocalIdx++,
- localScvIndices,
- true));
-
- globalScvfIndices_[eIdx].push_back(scvfIdx++);
- }
- }
- }
-
- // Compute the finite volume element geometries
- fvGeometries_.push_back(FVElementGeometry(*this, scvLocalIndexSet, scvfLocalIndexSet));
- }
- }
-
- // Binding of an element, has to be called before using the fvgeometries
- // Prepares all the volume variables within the element
- // For compatibility reasons with the FVGeometry cache being disabled
- void bindElement(const Element& element)
- {
- eIdxBound_ = problem_().elementMapper().index(element);
- }
- // this function is for compatibility reasons with cc methods
- void bind(const Element& element)
- {
- bindElement(element);
- }
-
-private:
- const Problem& problem_() const
- { return *problemPtr_; }
-
- const Problem* problemPtr_;
-
- IndexType eIdxBound_;
- GridView gridView_;
- Dumux::ElementMap<GridView> elementMap_;
- std::vector<std::shared_ptr<SubControlVolume>> scvs_;
- std::vector<std::shared_ptr<SubControlVolumeFace>> scvfs_;
- std::vector<FVElementGeometry> fvGeometries_;
-
- std::vector<std::vector<IndexType>> globalScvIndices_;
- std::vector<std::vector<IndexType>> globalScvfIndices_;
- const FeCache feCache_;
-};
-
-// specialization in case the FVElementGeometries are not stored
-template<class TypeTag>
-class BoxFVElementGeometryVector<TypeTag, false>
-{
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
- using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
-
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using CoordScalar = typename GridView::ctype;
-
- using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
- using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
- using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
-public:
- //! Constructor
- BoxFVElementGeometryVector(const GridView gridView)
- : gridView_(gridView), elementMap_(gridView), eIdxBound_(-1) {}
-
- /* \brief Get the finite volume geometry of an element
- * \note The finite volume geometry offers iterators over the sub control volumes
- * and the sub control volume faces of an element.
- */
- const FVElementGeometry& fvGeometry(const Element& element) const
- {
- auto eIdx = problem_().elementMapper().index(element);
- if (eIdx != eIdxBound_)
- DUNE_THROW(Dune::InvalidStateException, "Trying to get FVElementGeometry for element with index " << eIdx
- << ", but bound is element " << eIdxBound_ << ". Please call bindElement() before using the fvGeometry.");
- return fvGeometry_[0];
- }
-
- //! Get a sub control volume with a global scv index
- const SubControlVolume& subControlVolume(IndexType scvIdx) const
- {
- return stencilScvs_[scvIdx];
- }
-
- //! Get a sub control volume face with a global scvf index
- const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
- {
- return stencilScvfs_[scvfIdx];
- }
-
- //! Get a local finite element basis
- const FeLocalBasis& feLocalBasis(const Dune::GeometryType& type) const
- {
- return feCache_.get(type).localBasis();
- }
-
- //! The total number of sub control volumes
- std::size_t numScv() const
- {
- return numScvs_;
- }
-
- //! The total number of sun control volume faces
- std::size_t numScvf() const
- {
- return numScvf_;
- }
-
- //! The total number of boundary sub control volume faces
- //! For compatibility reasons with cc methods
- std::size_t numBoundaryScvf() const
- {
- return 0;
- }
-
- // Get an element from a sub control volume contained in it
- Element element(const SubControlVolume& scv) const
- { return elementMap_.element(scv.elementIndex()); }
-
- // Get an element from a global element index
- Element element(IndexType eIdx) const
- { return elementMap_.element(eIdx); }
-
- //! update all fvElementGeometries (do this again after grid adaption)
- void update(const Problem& problem)
- {
- problemPtr_ = &problem;
- elementMap_.clear();
- release_();
- eIdxBound_ = -1;
-
- auto numElems = gridView_.size(0);
- elementMap_.resize(numElems);
- scvIndices_.resize(numElems);
- scvFaceIndices_.resize(numElems);
-
- // save data on the grid's scvs and scvfs
- IndexType scvIdx = 0;
- IndexType scvfIdx = 0;
- for (const auto& element : elements(gridView_))
- {
- // the element-wise index sets for finite volume geometry
- std::vector<IndexType> scvIndexSet;
- std::vector<IndexType> scvfIndexSet;
-
- // fill the element map with seeds
- auto eIdx = problem.elementMapper().index(element);
- elementMap_[eIdx] = element.seed();
-
- // store scv indices of the element
- for (int scvLocalIdx = 0; scvLocalIdx < element.subEntities(dim); ++scvLocalIdx)
- scvIndexSet.push_back(scvIdx++);
-
- // store scv face indices of the element
- for (int scvfLocalIdx = 0; scvfLocalIdx < element.subEntities(dim-1); ++scvfLocalIdx)
- scvfIndexSet.push_back(scvfIdx++);
-
- // store the sub control volume face indices on the domain boundary
- for (const auto& intersection : intersections(gridView_, element))
- {
- if (intersection.boundary())
- {
- for (int localNodeIdx = 0; localNodeIdx < intersection.geometry().corners(); ++localNodeIdx)
- {
- // add this scvf to the element's scvf list
- scvfIndexSet.push_back(scvfIdx++);
- }
- }
- }
-
- // store the sets of indices in the data containers
- scvIndices_[eIdx] = scvIndexSet;
- scvFaceIndices_[eIdx] = scvfIndexSet;
- }
-
- numScvs_ = scvIdx;
- numScvf_ = scvfIdx;
- }
-
- // build the geometries inside an element
- void bindElement(const Element& element)
- {
- eIdxBound_ = problem_().elementMapper().index(element);
- release_();
- makeElementGeometries_(element);
- }
-
- // for cc methods, this method binds all the elements in the stencil.
- // Here, we simply forward to the bindElement() routine.
- void bind(const Element& element)
- {
- bindElement(element);
- }
-
-private:
- void release_()
- {
- stencilScvs_.clear();
- stencilScvfs_.clear();
- fvGeometry_.clear();
- }
-
- void makeElementGeometries_(const Element& element)
- {
- const auto eIdx = problem_().elementMapper().index(element);
-
- // the element-wise index sets for finite volume geometry
- auto scvIndexSet = scvIndices_[eIdx];
- auto scvfIndexSet = scvFaceIndices_[eIdx];
- std::vector<IndexType> localScvIndexSet;
- std::vector<IndexType> localScvfIndexSet;
- localScvIndexSet.reserve(scvIndexSet.size());
- localScvfIndexSet.reserve(scvfIndexSet.size());
-
- // get the element geometry
- auto elementGeometry = element.geometry();
- const auto& referenceElement = ReferenceElements::general(elementGeometry.type());
-
- // get the sub control volume geometries of this element
- auto boxGeometries = std::move(BoxGeometryHelper<GridView>::getScvAndScvfGeometries(elementGeometry));
- // define some aliases
- auto& scvGeometries = boxGeometries.first;
- auto& scvfGeometries = boxGeometries.second;
-
- // construct the sub control volumes
- IndexType scvLocalIdx = 0;
- for (auto&& scvGeometry : scvGeometries)
- {
- localScvIndexSet.push_back(scvLocalIdx);
- auto dofIdxGlobal = problem_().vertexMapper().subIndex(element, scvLocalIdx, dim);
-
- // add scv to the local container
- stencilScvs_.emplace_back(SubControlVolume(*scvGeometry,
- scvIndexSet[scvLocalIdx],
- eIdx,
- scvLocalIdx,
- dofIdxGlobal));
- // increment local counter
- scvLocalIdx++;
- }
-
- // construct the sub control volume faces
- IndexType scvfLocalIdx = 0;
- for (auto&& scvfGeometry : scvfGeometries)
- {
- localScvfIndexSet.push_back(scvfLocalIdx);
-
- // find the local scv indices this scvf is connsected to
- std::vector<IndexType> localScvIndices(2);
- localScvIndices[0] = referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim);
- localScvIndices[1] = referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim);
-
- // compute the scvf normal unit outer normal
- auto normal = std::move(BoxGeometryHelper<GridView>::normal(elementGeometry, *scvfGeometry));
- const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
- const auto s = v*normal;
- if (std::signbit(s))
- normal *= -1;
-
- stencilScvfs_.emplace_back(SubControlVolumeFace(*scvfGeometry,
- scvfGeometry->center(),
- normal,
- scvfIndexSet[scvfLocalIdx],
- scvfLocalIdx,
- localScvIndices,
- false));
-
- // increment local counter
- scvfLocalIdx++;
- }
-
- // construct the sub control volume faces on the domain boundary
- for (const auto& intersection : intersections(gridView_, element))
- {
- if (intersection.boundary())
- {
- auto isGeometry = intersection.geometry();
- auto boundaryScvfGeometries = std::move(BoxGeometryHelper<GridView>::getBoundaryScvfGeometries(isGeometry));
-
- IndexType boundaryFaceLocalIdx = 0;
- for (auto&& scvfGeometry : boundaryScvfGeometries)
- {
- localScvfIndexSet.push_back(scvfLocalIdx);
-
- // find the scv this scvf is connected to
- std::vector<IndexType> localScvIndices = {static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, boundaryFaceLocalIdx++, dim))};
-
- // get the unit outer normal through the intersection
- auto normal = intersection.unitOuterNormal(isGeometry.local(scvfGeometry->center()));
-
- stencilScvfs_.emplace_back(SubControlVolumeFace(*scvfGeometry,
- scvfGeometry->center(),
- normal,
- scvfIndexSet[scvfLocalIdx],
- scvfLocalIdx,
- localScvIndices,
- true));
-
- // increment local counter
- scvfLocalIdx++;
- }
- }
- }
-
- // Compute the finite volume element geometries
- fvGeometry_.push_back(FVElementGeometry(*this, localScvIndexSet, localScvfIndexSet));
- }
-
- const Problem& problem_() const
- { return *problemPtr_; }
-
- GridView gridView_;
- const Problem* problemPtr_;
- const FeCache feCache_;
-
- // Information on the global number of geometries
- IndexType numScvs_;
- IndexType numScvf_;
-
- // vectors that store the global data
- Dumux::ElementMap<GridView> elementMap_;
- std::vector<std::vector<IndexType>> scvIndices_;
- std::vector<std::vector<IndexType>> scvFaceIndices_;
-
- // vectors to store the geometries temporarily after binding an element
- IndexType eIdxBound_;
- std::vector<SubControlVolume> stencilScvs_;
- std::vector<SubControlVolumeFace> stencilScvfs_;
- std::vector<FVElementGeometry> fvGeometry_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/discretization/box/globalfvgeometry.hh b/dumux/discretization/box/globalfvgeometry.hh
new file mode 100644
index 0000000000..e1efe5803a
--- /dev/null
+++ b/dumux/discretization/box/globalfvgeometry.hh
@@ -0,0 +1,383 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the finite volume geometry vector for box models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element of the grid partition.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_GLOBAL_FVGEOMETRY_HH
+#define DUMUX_DISCRETIZATION_BOX_GLOBAL_FVGEOMETRY_HH
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dumux/discretization/box/boxgeometryhelper.hh>
+#include <dumux/discretization/box/fvelementgeometry.hh>
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/common/elementmap.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for box models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class BoxGlobalFVGeometry
+{};
+
+// specialization in case the FVElementGeometries are stored
+template<class TypeTag>
+class BoxGlobalFVGeometry<TypeTag, true>
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using Element = typename GridView::template Codim<0>::Entity;
+ //! The local class needs access to the scv, scvfs
+ //! as they are globally cached
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using CoordScalar = typename GridView::ctype;
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+
+ using GeometryHelper = BoxGeometryHelper<GridView, dim>;
+
+public:
+ //! Constructor
+ BoxGlobalFVGeometry(const GridView gridView)
+ : gridView_(gridView), elementMap_(gridView) {}
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ { return numScv_; }
+
+ //! The total number of sun control volume faces
+ std::size_t numScvf() const
+ { return numScvf_; }
+
+ //! The total number of boundary sub control volume faces
+ //! For compatibility reasons with cc methods
+ std::size_t numBoundaryScvf() const
+ { return numBoundaryScvf_; }
+
+ // Get an element from a sub control volume contained in it
+ Element element(const SubControlVolume& scv) const
+ { return elementMap_.element(scv.index()); }
+
+ // Get an element from a global element index
+ Element element(IndexType eIdx) const
+ { return elementMap_.element(eIdx); }
+
+ //! Return the gridView this global object lives on
+ const GridView& gridView() const
+ { return gridView_; }
+
+ //! update all fvElementGeometries (do this again after grid adaption)
+ void update(const Problem& problem)
+ {
+ problemPtr_ = &problem;
+
+ scvs_.clear();
+ scvfs_.clear();
+ elementMap_.clear();
+
+ auto numElements = gridView_.size(0);
+ scvs_.resize(numElements);
+ scvfs_.resize(numElements);
+ elementMap_.resize(numElements);
+
+ numScv_ = 0;
+ numScvf_ = 0;
+ numBoundaryScvf_ = 0;
+ // Build the SCV and SCV faces
+ for (const auto& element : elements(gridView_))
+ {
+ // fill the element map with seeds
+ auto eIdx = problem.elementMapper().index(element);
+ elementMap_[eIdx] = element.seed();
+
+ // count
+ numScv_ += element.subEntities(dim);
+ numScvf_ += element.subEntities(dim-1);
+
+ // get the element geometry
+ auto elementGeometry = element.geometry();
+ const auto& referenceElement = ReferenceElements::general(elementGeometry.type());
+
+ // get the sub control volume geometries of this element
+ GeometryHelper geometryHelper(elementGeometry);
+ auto scvGeometries = geometryHelper.createScvGeometries();
+ auto scvfGeometries = geometryHelper.createScvfGeometries();
+
+ // construct the sub control volumes
+ IndexType scvLocalIdx = 0;
+ scvs_[eIdx].reserve(scvGeometries.size());
+ for (auto&& scvGeometry : scvGeometries)
+ {
+ auto dofIdxGlobal = problem.vertexMapper().subIndex(element, scvLocalIdx, dim);
+ scvs_[eIdx].emplace_back(std::move(scvGeometry),
+ scvLocalIdx,
+ eIdx,
+ dofIdxGlobal);
+ // increment local counter
+ scvLocalIdx++;
+ }
+
+ // construct the sub control volume faces
+ IndexType scvfLocalIdx = 0;
+ scvfs_[eIdx].reserve(scvfGeometries.size());
+ for (auto&& scvfGeometry : scvfGeometries)
+ {
+ // find the global and local scv indices this scvf is belonging to
+ std::vector<IndexType> localScvIndices({static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+ static_cast<IndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+
+ // compute the scvf normal unit outer normal
+ auto normal = geometryHelper.normal(elementGeometry, *scvfGeometry);
+ const auto v = elementGeometry.corner(localScvIndices[1]) - elementGeometry.corner(localScvIndices[0]);
+ const auto s = v*normal;
+ if (std::signbit(s))
+ normal *= -1;
+
+ scvfs_.emplace_back(std::move(scvfGeometry),
+ normal,
+ scvfLocalIdx,
+ localScvIndices,
+ false);
+
+ // increment local counter
+ scvfLocalIdx++;
+ }
+
+ // construct the sub control volume faces on the domain boundary
+ for (const auto& intersection : intersections(gridView_, element))
+ {
+ if (intersection.boundary())
+ {
+ auto isGeometry = intersection.geometry();
+ // count
+ numScvf_ += isGeometry.corners();
+ numBoundaryScvf_ += isGeometry.corners();
+ auto boundaryScvfGeometries = geometryHelper.createBoundaryScvfGeometries(isGeometry);
+
+ IndexType isScvfLocalIdx = 0;
+ for (auto&& scvfGeometry : boundaryScvfGeometries)
+ {
+ // find the scvs this scvf is belonging to
+ std::vector<IndexType> localScvIndices =
+ {static_cast<IndexType>(referenceElement.subEntity(intersection.indexInInside(), 1,
+ isScvfLocalIdx++, dim))};
+
+ scvfs_.emplace_back(std::move(scvfGeometry),
+ intersection.centerUnitOuterNormal(),
+ scvfLocalIdx,
+ localScvIndices,
+ true);
+
+ // increment local counter
+ scvfLocalIdx++;
+ }
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend FVElementGeometry localView(const BoxGlobalFVGeometry& global)
+ {
+ return FVElementGeometry(global);
+ }
+
+ //! The finite element cache for creating local FE bases
+ const FeCache& feCache() const
+ { return feCache_; }
+
+private:
+
+ //! Get the local scvs for an element
+ const std::vector<SubControlVolume>& scvs(IndexType eIdx)
+ { return scvs_[eIdx]; }
+
+ //! Get the local scvfs for an element
+ const std::vector<SubControlVolume>& scvfs(IndexType eIdx)
+ { return scvfs_[eIdx]; }
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+ const FeCache feCache_;
+
+ GridView gridView_;
+ Dumux::ElementMap<GridView> elementMap_;
+ std::vector<std::vector<SubControlVolume>> scvs_;
+ std::vector<std::vector<SubControlVolumeFace>> scvfs_;
+ // TODO do we need those?
+ std::size_t numScv_;
+ std::size_t numScvf_;
+ std::size_t numBoundaryScvf_;
+};
+
+// specialization in case the FVElementGeometries are not stored
+template<class TypeTag>
+class BoxGlobalFVGeometry<TypeTag, false>
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ //! The local class needs access to the problem
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<dim>::Entity;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using CoordScalar = typename GridView::ctype;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
+
+public:
+ //! Constructor
+ BoxGlobalFVGeometry(const GridView gridView)
+ : gridView_(gridView), elementMap_(gridView) {}
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ { return numScv_; }
+
+ //! The total number of sun control volume faces
+ std::size_t numScvf() const
+ { return numScvf_; }
+
+ //! The total number of boundary sub control volume faces
+ //! For compatibility reasons with cc methods
+ std::size_t numBoundaryScvf() const
+ { return numBoundaryScvf_; }
+
+ // Get an element from a sub control volume contained in it
+ Element element(const SubControlVolume& scv) const
+ { return elementMap_.element(scv.index()); }
+
+ // Get an element from a global element index
+ Element element(IndexType eIdx) const
+ { return elementMap_.element(eIdx); }
+
+ //! Return the gridView this global object lives on
+ const GridView& gridView() const
+ { return gridView_; }
+
+ //! update all fvElementGeometries (do this again after grid adaption)
+ void update(const Problem& problem)
+ {
+ problemPtr_ = &problem;
+ elementMap_.clear();
+
+ auto numElems = gridView_.size(0);
+ elementMap_.resize(numElems);
+
+ // save global data on the grid's scvs and scvfs
+ numScv_ = 0;
+ numScvf_ = 0;
+ numBoundaryScvf_ = 0;
+ for (const auto& element : elements(gridView_))
+ {
+ // fill the element map with seeds
+ auto eIdx = problem.elementMapper().index(element);
+ elementMap_[eIdx] = element.seed();
+
+ numScv_ += element.subEntities(dim);
+ numScvf_ += element.subEntities(dim-1);
+
+ // store the sub control volume face indices on the domain boundary
+ for (const auto& intersection : intersections(gridView_, element))
+ {
+ if (intersection.boundary())
+ {
+ auto isGeometry = intersection.geometry();
+ numScvf_ += isGeometry.corners();
+ numBoundaryScvf_ += isGeometry.corners();
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend FVElementGeometry localView(const BoxGlobalFVGeometry& global)
+ {
+ return FVElementGeometry(global);
+ }
+
+ //! The finite element cache for creating local FE bases
+ const FeCache& feCache() const
+ { return feCache_; }
+
+private:
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ GridView gridView_;
+ const Problem* problemPtr_;
+ const FeCache feCache_;
+
+ // Information on the global number of geometries
+ // TODO do we need those information?
+ std::size_t numScv_;
+ std::size_t numScvf_;
+ std::size_t numBoundaryScvf_;
+
+ // vectors that store the global data
+ Dumux::ElementMap<GridView> elementMap_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/box/stencils.hh b/dumux/discretization/box/stencils.hh
index afe38dc3f8..489dead989 100644
--- a/dumux/discretization/box/stencils.hh
+++ b/dumux/discretization/box/stencils.hh
@@ -141,16 +141,16 @@ public:
elementStencils_[eIdx].update(problem, element);
// bind the FvGeometry to the element before using it
- problem.model().fvGeometries_().bindElement(element);
- const auto& fvGeometry = problem.model().fvGeometries(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
auto vIdxGlobal = scv.dofIndex();
vertexStencils_[vIdxGlobal].vertexScvs().push_back(scv.index());
vertexStencils_[vIdxGlobal].elementIndices().push_back(eIdx);
- for (const auto& scvJ : scvs(fvGeometry))
+ for (auto&& scvJ : scvs(fvGeometry))
{
auto vIdxGlobalJ = scvJ.dofIndex();
diff --git a/dumux/discretization/box/subcontrolvolume.hh b/dumux/discretization/box/subcontrolvolume.hh
index df886d74d3..410a2f272a 100644
--- a/dumux/discretization/box/subcontrolvolume.hh
+++ b/dumux/discretization/box/subcontrolvolume.hh
@@ -43,21 +43,15 @@ private:
public:
// the contructor in the box case
- BoxSubControlVolume(Geometry geometry,
- IndexType scvIdx,
- IndexType elementIndex,
- IndexType indexInElement,
- IndexType dofIndex)
+ BoxSubControlVolume(Geometry&& geometry,
+ IndexType scvIdx,
+ IndexType elementIndex,
+ IndexType dofIndex)
: SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)),
- scvIdx_(scvIdx), indexInElement_(std::move(indexInElement)),
+ scvIdx_(std::move(scvIdx)),
dofIndex_(std::move(dofIndex)) {}
- IndexType indexInElement() const
- {
- return indexInElement_;
- }
-
- //! The global index of this scv
+ //! The local index of this scv
IndexType index() const
{
return scvIdx_;
@@ -70,12 +64,11 @@ public:
GlobalPosition dofPosition() const
{
- return this->geometry().corner(indexInElement_);
+ return this->geometry().corner(scvIdx_);
}
private:
IndexType scvIdx_;
- IndexType indexInElement_;
IndexType dofIndex_;
};
diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh
index b96d5085c5..1f297a77e2 100644
--- a/dumux/discretization/box/subcontrolvolumeface.hh
+++ b/dumux/discretization/box/subcontrolvolumeface.hh
@@ -46,14 +46,12 @@ class BoxSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexT
using LocalPosition = Dune::FieldVector<Scalar, dim>;
public:
- BoxSubControlVolumeFace(const Geometry& geometry,
- const GlobalPosition& ipGlobal,
- const GlobalPosition& unitOuterNormal,
- IndexType scvfIndex,
- IndexType indexInElement,
- const std::vector<IndexType>& scvIndices,
- bool boundary = false)
- : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary) {}
+ BoxSubControlVolumeFace(Geometry&& geometry,
+ const GlobalPosition& unitOuterNormal,
+ IndexType scvfIndex,
+ const std::vector<IndexType>& scvIndices,
+ bool boundary = false)
+ : SubControlVolumeFaceBase<Geometry, IndexType>(std::move(geometry), geometry.center(), unitOuterNormal, scvfIndex, scvIndices, boundary) {}
};
} // end namespace
diff --git a/dumux/discretization/box/volumevariablesvector.hh b/dumux/discretization/box/volumevariablesvector.hh
index 313f39f796..31e88237bf 100644
--- a/dumux/discretization/box/volumevariablesvector.hh
+++ b/dumux/discretization/box/volumevariablesvector.hh
@@ -38,7 +38,7 @@ class BoxVolumeVariablesVector
// specialization in case of storing the volume variables
template<class TypeTag, bool useOldSol>
-class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true> : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables)>
+class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true>
{
friend BoxVolumeVariablesVector<TypeTag, !useOldSol, true>;
friend ImplicitModel<TypeTag>;
@@ -48,10 +48,10 @@ class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true> :
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
@@ -59,11 +59,8 @@ class BoxVolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true> :
BoxVolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const BoxVolumeVariablesVector<TypeTag, !useOldSol, true>& other)
{
- // do the copy
- numScvs_ = other.numScvs_;
- volumeVariables_ = other.volumeVariables_;
-
- // return the existing object
+ volVars_ = other.volVars_;
+ eIdx_ = other.eIdx_;
return *this;
};
@@ -72,52 +69,49 @@ public:
{
problemPtr_ = &problem;
- numScvs_ = problem.model().fvGeometries().numScv();
- volumeVariables_.resize(numScvs_);
+ volVars_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
- problem.model().fvGeometries_().bindElement(element);
- const auto& fvGeometry = problem.model().fvGeometries(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
(*this)[scv].update(sol[scv.dofIndex()],
- problem,
- element,
- scv);
+ problem,
+ element,
+ scv);
}
}
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
- const auto& scv = problem_().model().fvGeometries().subControlVolume(scvIdx);
- return (*this)[scv];
+ return volVars_[eIdx_][scvIdx];
}
VolumeVariables& operator [](IndexType scvIdx)
{
- const auto& scv = problem_().model().fvGeometries().subControlVolume(scvIdx);
- return (*this)[scv];
+ return volVars_[eIdx_][scvIdx];
}
const VolumeVariables& operator [](const SubControlVolume& scv) const
{
- return volumeVariables_[scv.index()];
+ return volVars_[eIdx_][scv.index()];
}
VolumeVariables& operator [](const SubControlVolume& scv)
{
- return volumeVariables_[scv.index()];
+ return volVars_[eIdx_][scv.index()];
}
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
- void bind(const Element& element) {}
- // In the box method, the vol vars within the elements adjacent to a vertex need to be bound
- void bind(const Vertex& vertex) {}
+ void bind(const Element& element)
+ { bindElement(element); }
// function to prepare the vol vars within the element
- void bindElement(const Element& element) {}
+ void bindElement(const Element& element)
+ { eIdx_ = problem_().elementMapper().index(element); }
private:
const Problem& problem_() const
@@ -125,9 +119,8 @@ private:
const Problem* problemPtr_;
- IndexType eIdxBound_;
- IndexType numScvs_;
- std::vector<VolumeVariables> volumeVariables_;
+ IndexType eIdx_;
+ std::vector<std::vector<VolumeVariables>> volVars_;
};
@@ -144,21 +137,16 @@ class BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
- BoxVolumeVariablesVector& operator= (const BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other) = default;
+ BoxVolumeVariablesVector& operator= (const BoxVolumeVariablesVector<TypeTag, false, false>& other) = default;
// operator curVolVars = prevVolVars
- void operator= (const BoxVolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>& other)
- {
- eIdxBound_ = -1;
- }
-
- BoxVolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
-
+ void operator= (const BoxVolumeVariablesVector<TypeTag, true, false>& other)
+ {}
public:
@@ -167,205 +155,54 @@ public:
problemPtr_ = &problem;
}
-
- // Binding of an element, prepares the volume variables within the element stencil
- // called by the local jacobian to prepare element assembly
- // specialization for cc models
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bind(const Element& element)
- {
- const auto eIdx = problem_().elementMapper().index(element);
- if (eIdx == eIdxBound_ && volVarIndices_.size() > 1)
- return;
-
- eIdxBound_ = eIdx;
-
- // make sure the FVElementGeometry is bound to the element
- problem_().model().fvGeometries_().bind(element);
- const auto& fvGeometry = problem_().model().fvGeometries(eIdx);
-
- // stencil information
- const auto& elementStencil = problem_().model().stencils(element).elementStencil();
- const auto& neighborStencil = problem_().model().stencils(element).neighborStencil();
-
- // maximum number of possible vol vars to be created
- const auto numDofs = elementStencil.size();;// + fvGeometry.numScvf();
-
- // resize local containers to the required size
- volumeVariables_.resize(numDofs);
- volVarIndices_.resize(numDofs);
- int localIdx = 0;
-
- // update the volume variables of the element at hand
- const auto& scvI = problem_().model().fvGeometries().subControlVolume(eIdx);
- const auto& solI = problem_().model().curSol()[eIdx];
- // VolumeVariables tmp;
- // tmp.update(solI, problem_(), element, scvI);
- volumeVariables_[localIdx].update(solI, problem_(), element, scvI);
- volVarIndices_[localIdx] = scvI.index();
- // volumeVariables_.push_back(tmp);
- // volVarIndices_.push_back(scvI.index());
- localIdx++;
-
- // Update the volume variables of the neighboring elements
- for (auto globalJ : neighborStencil)
- {
- const auto& elementJ = problem_().model().fvGeometries().element(globalJ);
- const auto& scvJ = problem_().model().fvGeometries().subControlVolume(globalJ);
- const auto& solJ = problem_().model().curSol()[globalJ];
- // tmp.update(solJ, problem_(), elementJ, scvJ);
- // volumeVariables_.push_back(tmp);
- // volVarIndices_.push_back(scvJ.index());
- volumeVariables_[localIdx].update(solJ, problem_(), elementJ, scvJ);
- volVarIndices_[localIdx] = scvJ.index();
- localIdx++;
- }
-
- // Update boundary volume variables
- for (const auto& scvFace : scvfs(fvGeometry))
- {
- // if we are not on a boundary, skip the rest
- if (!scvFace.boundary())
- continue;
-
- // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries
- const auto bcTypes = problem_().boundaryTypes(element, scvFace);
- if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */bcTypes.hasNeumann() || bcTypes.hasOutflow())
- continue;
-
- const auto dirichletPriVars = problem_().dirichlet(element, scvFace);
- // tmp.update(dirichletPriVars, problem_(), element, scvI);
- // volumeVariables_.push_back(tmp);
- // volVarIndices_.push_back(scvI.index());
- volumeVariables_.resize(localIdx+1);
- volVarIndices_.resize(localIdx+1);
- volumeVariables_[localIdx].update(dirichletPriVars, problem_(), element, scvI);
- volVarIndices_[localIdx] = scvFace.outsideScvIdx();
- localIdx++;
- }
- }
-
// specialization for box models, simply forwards to the bindElement method
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bind(const Element& element)
+ void bind(const Element& element, const FVElementGeometry& fvGeometry)
{
- bindElement(element);
+ bindElement(element, fvGeometry);
}
- // Binding of an element, prepares only the volume variables of the element
- // specialization for cc models
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bindElement(const Element& element)
- {
- const auto eIdx = problem_().elementMapper().index(element);
- if (eIdx == eIdxBound_ || std::find(volVarIndices_.begin(), volVarIndices_.end(), eIdx) != volVarIndices_.end())
- return;
-
- volumeVariables_.resize(1);
- volVarIndices_.resize(1);
- eIdxBound_ = eIdx;
-
- // make sure the FVElementGeometry is bound to the element
- problem_().model().fvGeometries_().bindElement(element);
-
- // update the volume variables of the element at hand
- const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
- const auto& sol = problem_().model().curSol()[eIdx];
- volumeVariables_[0].update(sol, problem_(), element, scv);
- volVarIndices_[0] = scv.index();
- }
-
- // In the box method, the vol vars within the elements adjacent to a vertex need to be bound (TODO: IMPLEMENT)
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bind(const Vertex& vertex) const {}
-
// specialization for box models
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bindElement(const Element& element)
+ void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
{
- const auto eIdx = problem_().elementMapper().index(element);
- if (eIdx == eIdxBound_)
- return;
-
- eIdxBound_ = eIdx;
-
- // make sure the FVElementGeometry is bound to the element
- problem_().model().fvGeometries_().bind(element);
- const auto& fvGeometry = problem_().model().fvGeometries(element);
-
- // get stencil information
- const auto numDofs = element.subEntities(dim);
-
// resize volume variables to the required size
- volumeVariables_.resize(numDofs);
- // volVarIndices_.resize(numDofs);
+ volVars_.resize(fvGeometry.numScv());
- int localIdx = 0;
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
- // std::cout << "scv index: " << scv.index() << ", dofIdx: " << scv.dofIndex() << ", localIdx: " << scv.indexInElement() << std::endl;
const auto& sol = problem_().model().curSol()[scv.dofIndex()];
- // let the interface solver update the volvars
- // volVarIndices_[localIdx] = scv.index();
// TODO: INTERFACE SOLVER
// problem_().model().boxInterfaceConditionSolver().updateScvVolVars(element, scv, sol);
- volumeVariables_[scv.indexInElement()].update(sol, problem_(), element, scv);
- localIdx++;
+ volVars_[scv.index()].update(sol, problem_(), element, scv);
}
- // std::cout << "finished\n";
}
const VolumeVariables& operator [](IndexType scvIdx) const
{
- return volumeVariables_[scvIdx];
+ return volVars_[scvIdx];
}
VolumeVariables& operator [](IndexType scvIdx)
{
- return volumeVariables_[scvIdx];
+ return volVars_[scvIdx];
}
const VolumeVariables& operator [](const SubControlVolume& scv) const
{
- return volumeVariables_[scv.indexInElement()];
+ return volVars_[scv.index()];
}
VolumeVariables& operator [](const SubControlVolume& scv)
{
- return volumeVariables_[scv.indexInElement()];
+ return volVars_[scv.index()];
}
private:
- void release_()
- {
- volumeVariables_.clear();
- volVarIndices_.clear();
- }
-
- const int getLocalIdx_(const int volVarIdx) const
- {
- auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
-
- if (it != volVarIndices_.end())
- return std::distance(volVarIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not find the current volume variables for volVarIdx = " << volVarIdx <<
- ", make sure to properly bind the volume variables to the element before using them");
- }
-
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
- IndexType eIdxBound_;
- std::vector<IndexType> volVarIndices_;
- std::vector<VolumeVariables> volumeVariables_;
+ std::vector<VolumeVariables> volVars_;
};
// Specialization when the previous volume variables are not stored
@@ -381,22 +218,18 @@ class BoxVolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using IndexType = typename GridView::IndexSet::IndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
BoxVolumeVariablesVector& operator= (const BoxVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other)
{
- release_();
+ volVars_.clear();
problemPtr_ = other.problemPtr_;
- eIdxBound_ = -1;
return *this;
}
- BoxVolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
-
-
public:
void update(const Problem& problem, const SolutionVector& sol)
@@ -404,116 +237,49 @@ public:
problemPtr_ = &problem;
}
- // Binding of an element, prepares the volume variables of only the element
- // specialization for cc models
- template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bindElement(const Element& element)
- {
- const auto eIdx = problem_().elementMapper().index(element);
- if (eIdx == eIdxBound_)
- return;
-
- volumeVariables_.resize(1);
- volVarIndices_.resize(1);
- eIdxBound_ = eIdx;
-
- // make sure FVElementGeometry is bound to the element
- problem_().model().fvGeometries_().bindElement(element);
-
- // update the volume variables of the element at hand
- const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
- const auto& sol = problem_().model().prevSol()[eIdx];
- volumeVariables_[0].update(sol, problem_(), element, scv);
- volVarIndices_[0] = scv.index();
- }
-
- // specialization for box models
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bindElement(const Element& element)
+ // bind an element
+ void bindElement(const Element& element, const FVElementGeometry& fvGeometry)
{
- const auto eIdx = problem_().elementMapper().index(element);
- if (eIdx == eIdxBound_)
- return;
-
- eIdxBound_ = eIdx;
-
- // make sure the FVElementGeometry is bound to the element
- problem_().model().fvGeometries_().bind(element);
- const auto& fvGeometry = problem_().model().fvGeometries(element);
-
- // get stencil information
- const auto numDofs = element.subEntities(dim);
-
// resize volume variables to the required size
- volumeVariables_.resize(numDofs);
- // volVarIndices_.resize(numDofs);
+ volVars_.resize(fvGeometry.numScv());
- int localIdx = 0;
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
const auto& sol = problem_().model().prevSol()[scv.dofIndex()];
- // let the interface solver update the volvars
- // volVarIndices_[localIdx] = scv.index();
//TODO: INTERFACE SOLVER?
// problem_().model().boxInterfaceConditionSolver().updateScvVolVars(element, scv, sol);
- volumeVariables_[localIdx].update(sol, problem_(), element, scv);
- localIdx++;
+ volVars_[scv.index()].update(sol, problem_(), element, scv);
}
}
- // In the box method, the vol vars within the elements adjacent to a vertex need to be bound (TODO: IMPLEMENT)
- template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox)>::type
- bind(const Vertex& vertex) const {}
-
- // const VolumeVariables& operator [](IndexType scvIdx) const
- // {
- // return volumeVariables_[getLocalIdx_(scvIdx)];
- // }
-
- // VolumeVariables& operator [](IndexType scvIdx)
- // {
- // return volumeVariables_[getLocalIdx_(scvIdx)];
- // }
-
- const VolumeVariables& operator [](const SubControlVolume& scv) const
+ const VolumeVariables& operator [](IndexType scvIdx) const
{
- return volumeVariables_[scv.indexInElement()];
+ return volVars_[scvIdx];
}
- VolumeVariables& operator [](const SubControlVolume& scv)
+ VolumeVariables& operator [](IndexType scvIdx)
{
- return volumeVariables_[scv.indexInElement()];
+ return volVars_[scvIdx];
}
-private:
-
- void release_()
+ const VolumeVariables& operator [](const SubControlVolume& scv) const
{
- volumeVariables_.clear();
- volVarIndices_.clear();
+ return volVars_[scv.index()];
}
- const int getLocalIdx_(const int volVarIdx) const
+ VolumeVariables& operator [](const SubControlVolume& scv)
{
- auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
-
- if (it != volVarIndices_.end())
- return std::distance(volVarIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not find the previous volume variables for volVarIdx = " << volVarIdx <<
- ", make sure to properly bind the volume variables to the element before using them");
+ return volVars_[scv.index()];
}
+private:
+
Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
- IndexType eIdxBound_;
- std::vector<IndexType> volVarIndices_;
- std::vector<VolumeVariables> volumeVariables_;
+
+ std::vector<VolumeVariables> volVars_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/fluxvariablescachevector.hh b/dumux/discretization/cellcentered/fluxvariablescachevector.hh
index 859f16d12f..a44c8084f4 100644
--- a/dumux/discretization/cellcentered/fluxvariablescachevector.hh
+++ b/dumux/discretization/cellcentered/fluxvariablescachevector.hh
@@ -39,24 +39,26 @@ class CCFluxVariablesCacheVector
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
// When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ update(Problem& problem)
{
problemPtr_ = &problem;
-
- fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
+ const auto& globalFvGeometry = problem.model().globalFvGeometry();
+ fluxVarsCache_.resize(globalFvGeometry.numScvf());
for (const auto& element : elements(problem.gridView()))
{
// Prepare the geometries within the elements of the stencil
- problem.model().fvGeometries_().bind(element);
+ auto fvGeometry = localView(globalFvGeometry);
+ fvGeometry.bind(element);
- const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
(*this)[scvf].update(problem, element, scvf);
}
@@ -65,7 +67,8 @@ public:
// When global flux variables caching is disabled, we don't need to update the cache
template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type update(Problem& problem)
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ update(Problem& problem)
{
problemPtr_ = &problem;
}
@@ -73,10 +76,9 @@ public:
// This function has to be called prior to flux calculations on the element.
// Prepares the transmissibilities of the scv faces in an element. The FvGeometry is assumed to be bound.
template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bindElement(const Element& element)
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ bindElement(const Element& element, const FVElementGeometry& fvGeometry)
{
- const auto& fvGeometry = problem_().model().fvGeometries(element);
-
// resizing of the cache
const auto numScvf = fvGeometry.numScvf();
fluxVarsCache_.resize(numScvf);
@@ -84,9 +86,9 @@ public:
IndexType localScvfIdx = 0;
// fill the containers
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
- fluxVarsCache_[localScvfIdx].update(problem_(), element, scvf);
+ fluxVarsCache_[localScvfIdx].update(problem_(), element, fvGeometry, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
@@ -94,20 +96,22 @@ public:
// Specialization for the global caching being enabled - do nothing here
template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bindElement(const Element& element) {}
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ bindElement(const Element& element, const FVElementGeometry& fvGeometry)
+ {}
// This function is called by the CCLocalResidual before flux calculations during assembly.
// Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound.
template <typename T = TypeTag>
- typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bind(const Element& element)
+ typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ bind(const Element& element, const FVElementGeometry& fvGeometry)
{
const auto globalI = problem_().elementMapper().index(element);
- const auto& fvGeometry = problem_().model().fvGeometries(element);
const auto& neighborStencil = problem_().model().stencils(element).neighborStencil();
const auto numNeighbors = neighborStencil.size();
// find the number of scv faces that need to be prepared
- IndexType numScvf = fvGeometry.numScvf();
+ auto numScvf = fvGeometry.numScvf();
for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
{
const auto& fluxVarIndicesJ = problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
@@ -118,9 +122,9 @@ public:
fluxVarsCache_.resize(numScvf);
globalScvfIndices_.resize(numScvf);
IndexType localScvfIdx = 0;
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
- fluxVarsCache_[localScvfIdx].update(problem_(), element, scvf);
+ fluxVarsCache_[localScvfIdx].update(problem_(), element, fvGeometry, scvf);
globalScvfIndices_[localScvfIdx] = scvf.index();
localScvfIdx++;
}
@@ -129,11 +133,12 @@ public:
for (IndexType localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
{
const auto& fluxVarIndicesJ = problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
- const auto elementJ = problem_().model().fvGeometries().element(neighborStencil[localIdxJ]);
+ const auto elementJ = fvGeometry.globalFvGeometry().element(neighborStencil[localIdxJ]);
+
for (auto fluxVarIdx : fluxVarIndicesJ)
{
- const auto& scvfJ = problem_().model().fvGeometries().subControlVolumeFace(fluxVarIdx);
- fluxVarsCache_[localScvfIdx].update(problem_(), elementJ, scvfJ);
+ auto&& scvfJ = fvGeometry.scvf(fluxVarIdx);
+ fluxVarsCache_[localScvfIdx].update(problem_(), elementJ, fvGeometry, scvfJ);
globalScvfIndices_[localScvfIdx] = scvfJ.index();
localScvfIdx++;
}
@@ -142,7 +147,8 @@ public:
// Specialization for the global caching being enabled - do nothing here
template <typename T = TypeTag>
- typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type bind(const Element& element) {}
+ typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalFluxVariablesCache)>::type
+ bind(const Element& element, const FVElementGeometry& fvGeometry) {}
// access operators in the case of caching
template <typename T = TypeTag>
@@ -173,11 +179,8 @@ private:
getLocalScvfIdx_(const int scvfIdx) const
{
auto it = std::find(globalScvfIndices_.begin(), globalScvfIndices_.end(), scvfIdx);
-
- if (it != globalScvfIndices_.end())
- return std::distance(globalScvfIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not find the flux vars cache for scvfIdx = " << scvfIdx);
+ assert(it != globalScvfIndices_.end() && "Could not find the flux vars cache for scvfIdx");
+ return std::distance(globalScvfIndices_.begin(), it);
}
const Problem& problem_() const
diff --git a/dumux/discretization/cellcentered/stencils.hh b/dumux/discretization/cellcentered/stencils.hh
index 98c0a346d8..ef2b34e86c 100644
--- a/dumux/discretization/cellcentered/stencils.hh
+++ b/dumux/discretization/cellcentered/stencils.hh
@@ -38,22 +38,24 @@ class CCElementStencils
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
using Element = typename GridView::template Codim<0>::Entity;
// TODO a separate stencil class all stencils can derive from?
using Stencil = std::vector<IndexType>;
public:
- void update(const Problem& problem, const Element& element)
+ //! Update the stencil. We expect a bound fvGeometry
+ void update(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry)
{
elementStencil_.clear();
- const auto& fvGeometry = problem.model().fvGeometries(element);
// loop over sub control faces
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
FluxVariables fluxVars;
- const auto& stencil = fluxVars.computeStencil(problem, element, scvf);
-
+ const auto& stencil = fluxVars.computeStencil(problem, element, fvGeometry, scvf);
elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end());
}
// make values in elementstencil unique
@@ -64,15 +66,16 @@ public:
neighborStencil_ = elementStencil_;
// remove the element itself and possible ghost neighbors from the neighbor stencil
- auto pred = [&problem, globalI](const int i) -> bool
+ auto pred = [&fvGeometry, globalI](const int i) -> bool
{
if (i == globalI)
return true;
- if (problem.model().fvGeometries().element(i).partitionType() == Dune::GhostEntity)
+ if (fvGeometry.globalFvGeometry().element(i).partitionType() == Dune::GhostEntity)
return true;
return false;
};
- neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), pred), neighborStencil_.end());
+ neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), pred),
+ neighborStencil_.end());
}
//! The full element stencil (all element this element is interacting with)
@@ -110,11 +113,12 @@ public:
elementStencils_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
- // bind the FvGeometry to the element before using it
- problem.model().fvGeometries_().bind(element);
+ // restrict the FvGeometry locally and bind to the element
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bind(element);
auto eIdx = problem.elementMapper().index(element);
- elementStencils_[eIdx].update(problem, element);
+ elementStencils_[eIdx].update(problem, element, fvGeometry);
}
}
diff --git a/dumux/discretization/cellcentered/tpfa/CMakeLists.txt b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt
index ff7ebeb2ba..9bb40dd863 100644
--- a/dumux/discretization/cellcentered/tpfa/CMakeLists.txt
+++ b/dumux/discretization/cellcentered/tpfa/CMakeLists.txt
@@ -1,5 +1,6 @@
install(FILES
-fvelementgeometryvector.hh
+globalfvgeometry.hh
+fvelementgeometry.hh
subcontrolvolumeface.hh
subcontrolvolume.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/discretization/cellcentered/tpfa)
\ No newline at end of file
diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index c24d224618..345e94d58c 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -52,33 +52,35 @@ NEW_PROP_TAG(ProblemEnableGravity);
template <class TypeTag>
class DarcysLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, DiscretizationMethod) == GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa>::type >
{
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::IndexSet::IndexType IndexType;
- typedef std::vector<IndexType> Stencil;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Stencil = std::vector<IndexType>;
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
- typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimWorldMatrix;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
static Scalar flux(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvFace,
const IndexType phaseIdx)
{
const auto& tij = problem.model().fluxVarsCache(scvFace).tij();
// Get the inside and outside volume variables
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
+ const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
const auto& insideVolVars = problem.model().curVolVars(insideScv);
const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
@@ -108,7 +110,9 @@ public:
else
{
const auto outsideScvIdx = scvFace.outsideScvIdx();
- const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+ // as we assemble fluxes from the neighbor to our element the outside index
+ // refers to the scv of our element, so we use the scv method
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto xOutside = outsideScv.center();
const auto gOutside = problem.gravityAtPos(xOutside);
hOutside -= rho*(gOutside*xOutside);
@@ -118,7 +122,10 @@ public:
return tij*(hInside - hOutside);
}
- static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+ static Stencil stencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
{
Stencil stencil;
if (!scvFace.boundary())
@@ -134,20 +141,25 @@ public:
// The flux variables cache has to be bound to an element prior to flux calculations
// During the binding, the transmissibilities will be computed and stored using the method below.
- static Scalar calculateTransmissibilities(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+ static Scalar calculateTransmissibilities(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
{
Scalar tij;
const auto insideScvIdx = scvFace.insideScvIdx();
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv);
Scalar ti = calculateOmega_(problem, scvFace, insideK, element, insideScv);
if (!scvFace.boundary())
{
const auto outsideScvIdx = scvFace.outsideScvIdx();
- const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
- const auto outsideElement = problem.model().fvGeometries().element(outsideScvIdx);
+ // as we assemble fluxes from the neighbor to our element the outside index
+ // refers to the scv of our element, so we use the scv method
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
+ const auto outsideElement = fvGeometry.globalFvGeometry().element(outsideScvIdx);
const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv);
Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideElement, outsideScv);
@@ -163,7 +175,11 @@ public:
private:
- static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &K, const Element& element, const SubControlVolume &scv)
+ static Scalar calculateOmega_(const Problem& problem,
+ const SubControlVolumeFace& scvFace,
+ const DimWorldMatrix &K,
+ const Element& element,
+ const SubControlVolume &scv)
{
GlobalPosition Knormal;
K.mv(scvFace.unitOuterNormal(), Knormal);
@@ -178,7 +194,11 @@ private:
return omega;
}
- static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const Scalar K, const Element& element, const SubControlVolume &scv)
+ static Scalar calculateOmega_(const Problem& problem,
+ const SubControlVolumeFace& scvFace,
+ const Scalar K,
+ const Element& element,
+ const SubControlVolume &scv)
{
auto distanceVector = scvFace.center();
distanceVector -= scv.center();
diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
index cddca62731..6f6d5cc37b 100644
--- a/dumux/discretization/cellcentered/tpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
@@ -59,6 +59,7 @@ class FicksLaw<TypeTag, typename std::enable_if<GET_PROP_VALUE(TypeTag, Discreti
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::IndexSet::IndexType IndexType;
typedef typename std::vector<IndexType> Stencil;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
@@ -73,16 +74,16 @@ public:
static Scalar flux(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvFace,
const int phaseIdx,
const int compIdx)
{
// diffusion tensors are always solution dependent
- Scalar tij = calculateTransmissibility_(problem, element, scvFace, phaseIdx, compIdx);
+ Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, scvFace, phaseIdx, compIdx);
// Get the inside volume variables
- const auto insideScvIdx = scvFace.insideScvIdx();
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+ const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
const auto& insideVolVars = problem.model().curVolVars(insideScv);
// and the outside volume variables
@@ -96,7 +97,10 @@ public:
return rho*tij*(xInside - xOutside);
}
- static Stencil stencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+ static Stencil stencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
{
std::vector<IndexType> stencil;
stencil.clear();
@@ -114,12 +118,16 @@ public:
private:
- static Scalar calculateTransmissibility_(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx)
+ static Scalar calculateTransmissibility_(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace,
+ const int phaseIdx, const int compIdx)
{
Scalar tij;
const auto insideScvIdx = scvFace.insideScvIdx();
- const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
@@ -129,7 +137,7 @@ private:
if (!scvFace.boundary())
{
const auto outsideScvIdx = scvFace.outsideScvIdx();
- const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+ const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
new file mode 100644
index 0000000000..9acef62f15
--- /dev/null
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -0,0 +1,360 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for a local finite volume geometry for cell-centered TPFA models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element in the local scope we are restricting to, e.g. stencil or element.
+ */
+#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CCTPFA_FV_ELEMENT_GEOMETRY_HH
+
+#include <dune/common/iteratorrange.hh>
+
+#include <dumux/discretization/scvandscvfiterators.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+
+namespace Dumux
+{
+
+//! forward declaration of the global finite volume geometry
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class CCTpfaGlobalFVGeometry;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for cell-centered TPFA models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class CCTpfaFVElementGeometry
+{};
+
+//! specialization in case the FVElementGeometries are stored globally
+//! In this case we just forward internally to the global object
+template<class TypeTag>
+class CCTpfaFVElementGeometry<TypeTag, true>
+{
+ using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+
+ using ScvIterator = ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
+ using ScvfIterator = ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
+
+public:
+ //! Constructor
+ CCTpfaFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
+ : globalFvGeometryPtr_(&globalFvGeometry) {}
+
+ //! Get an elment sub control volume with a global scv index
+ //! We separate element and neighbor scvs to speed up mapping
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ return globalFvGeometry().scv(scvIdx);
+ }
+
+ //! Get an element sub control volume face with a global scvf index
+ //! We separate element and neighbor scvfs to speed up mapping
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ return globalFvGeometry().scvf(scvfIdx);
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element (not including neighbor scvs)
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline Dune::IteratorRange<ScvIterator>
+ scvs(const CCTpfaFVElementGeometry& fvGeometry)
+ {
+ return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
+ ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element (not including neighbor scvfs)
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline Dune::IteratorRange<ScvfIterator>
+ scvfs(const CCTpfaFVElementGeometry& fvGeometry)
+ {
+ const auto& g = fvGeometry.globalFvGeometry();
+ const auto scvIdx = fvGeometry.scvIndices()[0];
+ return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry),
+ ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry));
+ }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScv() const
+ {
+ return scvIndices_.size();
+ }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScvf() const
+ {
+ return globalFvGeometry().scvfIndicesOfScv(scvIndices_[0]).size();
+ }
+
+ //! Binding of an element, called by the local jacobian to prepare element assembly
+ void bind(const Element& element)
+ {
+ this->bindElement(element);
+ }
+
+ //! Bind only element-local
+ void bindElement(const Element& element)
+ {
+ elementPtr_ = &element;
+ scvIndices_ = std::vector<IndexType>({globalFvGeometry().problem_().elementMapper().index(*elementPtr_)});
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ const GlobalFVGeometry& globalFvGeometry() const
+ { return *globalFvGeometryPtr_; }
+
+private:
+
+ const Element* elementPtr_;
+ std::vector<IndexType> scvIndices_;
+ const GlobalFVGeometry* globalFvGeometryPtr_;
+};
+
+//! specialization in case the FVElementGeometries are not stored
+template<class TypeTag>
+class CCTpfaFVElementGeometry<TypeTag, false>
+{
+ using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+
+ using ScvIterator = ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
+ using ScvfIterator = ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
+
+public:
+ //! Constructor
+ CCTpfaFVElementGeometry(const GlobalFVGeometry& globalFvGeometry)
+ : globalFvGeometryPtr_(&globalFvGeometry) {}
+
+ //! Get an elment sub control volume with a global scv index
+ //! We separate element and neighbor scvs to speed up mapping
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ if (scvIdx == scvIndices_[0])
+ return scvs_[0];
+ else
+ return neighborScvs_[findLocalIndex(scvIdx, neighborScvIndices_)];
+ }
+
+ //! Get an element sub control volume face with a global scvf index
+ //! We separate element and neighbor scvfs to speed up mapping
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
+ if (it != scvfIndices_.end())
+ return scvfs_[std::distance(scvfIndices_.begin(), it)];
+ else
+ return neighborScvfs_[findLocalIndex(scvfIdx, neighborScvfIndices_)];
+ }
+
+ //! iterator range for sub control volumes. Iterates over
+ //! all scvs of the bound element (not including neighbor scvs)
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volumes of this FVElementGeometry use
+ //! for (auto&& scv : scvs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
+ scvs(const ThisType& g)
+ {
+ using Iter = typename std::vector<SubControlVolume>::const_iterator;
+ return Dune::IteratorRange<Iter>(g.scvs_.begin(), g.scvs_.end());
+ }
+
+ //! iterator range for sub control volumes faces. Iterates over
+ //! all scvfs of the bound element (not including neighbor scvfs)
+ //! This is a free function found by means of ADL
+ //! To iterate over all sub control volume faces of this FVElementGeometry use
+ //! for (auto&& scvf : scvfs(fvGeometry))
+ friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
+ scvfs(const ThisType& g)
+ {
+ using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<Iter>(g.scvfs_.begin(), g.scvfs_.end());
+ }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScv() const
+ {
+ return scvs_.size();
+ }
+
+ //! number of sub control volumes in this fv element geometry
+ std::size_t numScvf() const
+ {
+ return scvfs_.size();
+ }
+
+ //! Binding of an element preparing the geometries of the whole stencil
+ //! called by the local jacobian to prepare element assembly
+ void bind(const Element& element)
+ {
+ bindElement(element);
+ for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
+ {
+ if (intersection.neighbor())
+ makeNeighborGeometries(intersection.outside());
+ }
+ }
+
+ //! Binding of an element preparing the geometries only inside the element
+ void bindElement(const Element& element)
+ {
+ clear();
+ elementPtr_ = &element;
+ makeElementGeometries(element);
+ }
+
+ //! The global finite volume geometry we are a restriction of
+ const GlobalFVGeometry& globalFvGeometry() const
+ { return *globalFvGeometryPtr_; }
+
+private:
+
+ //! create scvs and scvfs of the bound element
+ void makeElementGeometries(const Element& element)
+ {
+ auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
+ scvs_.emplace_back(element.geometry(), eIdx);
+ scvIndices_.emplace_back(eIdx);
+
+ const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx);
+ const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx);
+
+ int scvfCounter = 0;
+ for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
+ {
+ if (intersection.neighbor() || intersection.boundary())
+ {
+ scvfs_.emplace_back(intersection,
+ scvFaceIndices[scvfCounter],
+ std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
+ );
+
+ scvfIndices_.push_back(scvFaceIndices[scvfCounter]);
+ scvfCounter++;
+ }
+ }
+ }
+
+ //! create the necessary scvs and scvfs of the neighbor elements to the bound elements
+ void makeNeighborGeometries(const Element& element)
+ {
+ // create the neighbor scv
+ auto eIdx = globalFvGeometry().problem_().elementMapper().index(element);
+ neighborScvs_.emplace_back(element.geometry(), eIdx);
+ neighborScvIndices_.push_back(eIdx);
+
+ const auto& scvFaceIndices = globalFvGeometry().scvfIndicesOfScv(eIdx);
+ const auto& neighborVolVarIndices = globalFvGeometry().neighborVolVarIndices(eIdx);
+
+ int scvfCounter = 0;
+ for (const auto& intersection : intersections(globalFvGeometry().gridView(), element))
+ {
+ if (intersection.neighbor())
+ {
+ // only create subcontrol faces where the outside element is the bound element
+ if (intersection.outside() == *elementPtr_)
+ {
+ neighborScvfs_.emplace_back(intersection,
+ scvFaceIndices[scvfCounter],
+ std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
+ );
+
+ neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
+ }
+ // always increase local counter
+ scvfCounter++;
+ }
+ else if (intersection.boundary())
+ {
+ neighborScvfs_.emplace_back(intersection,
+ scvFaceIndices[scvfCounter],
+ std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
+ );
+
+ neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
+ scvfCounter++;
+ }
+ }
+ }
+
+ const IndexType findLocalIndex(const IndexType idx,
+ const std::vector<IndexType>& indices) const
+ {
+ auto it = std::find(indices.begin(), indices.end(), idx);
+ assert(it != indices.end() && "Could not find the scv/scvf! Make sure to properly bind this class!");
+ return std::distance(indices.begin(), it);
+
+ }
+
+ void clear()
+ {
+ scvIndices_.clear();
+ scvfIndices_.clear();
+ scvs_.clear();
+ scvfs_.clear();
+
+ neighborScvIndices_.clear();
+ neighborScvfIndices_.clear();
+ neighborScvs_.clear();
+ neighborScvfs_.clear();
+ }
+
+ // the bound element
+ const Element* elementPtr_;
+
+ const GlobalFVGeometry* globalFvGeometryPtr_;
+
+ // local storage after binding an element
+ std::vector<IndexType> scvIndices_;
+ std::vector<IndexType> scvfIndices_;
+ std::vector<SubControlVolume> scvs_;
+ std::vector<SubControlVolumeFace> scvfs_;
+
+ std::vector<IndexType> neighborScvIndices_;
+ std::vector<IndexType> neighborScvfIndices_;
+ std::vector<SubControlVolume> neighborScvs_;
+ std::vector<SubControlVolumeFace> neighborScvfs_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh
similarity index 59%
rename from dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh
rename to dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh
index 1a239b4e38..35e3e1f15d 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh
+++ b/dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh
@@ -20,13 +20,14 @@
* \file
* \brief Base class for the finite volume geometry vector for cell-centered TPFA models
* This builds up the sub control volumes and sub control volume faces
- * for each element.
+ * for each element of the grid partition.
*/
-#ifndef DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH
-#define DUMUX_DISCRETIZATION_CCTPFA_FV_GEOMETRY_VECTOR_HH
+#ifndef DUMUX_DISCRETIZATION_CCTPFA_GLOBAL_FVGEOMETRY_HH
+#define DUMUX_DISCRETIZATION_CCTPFA_GLOBAL_FVGEOMETRY_HH
-#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/common/elementmap.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
namespace Dumux
{
@@ -37,13 +38,13 @@ namespace Dumux
* This builds up the sub control volumes and sub control volume faces
* for each element.
*/
-template<class TypeTag, bool EnableFVElementGeometryCache>
-class CCTpfaFVElementGeometryVector
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class CCTpfaGlobalFVGeometry
{};
-// specialization in case the FVElementGeometries are stored
+// specialization in case the FVElementGeometries are stored globally
template<class TypeTag>
-class CCTpfaFVElementGeometryVector<TypeTag, true>
+class CCTpfaGlobalFVGeometry<TypeTag, true>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -52,33 +53,15 @@ class CCTpfaFVElementGeometryVector<TypeTag, true>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
+ //! The local class needs access to the scv, scvfs and the fv element geometry
+ //! as they are globally cached
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
public:
//! Constructor
- CCTpfaFVElementGeometryVector(const GridView gridView)
+ CCTpfaGlobalFVGeometry(const GridView gridView)
: gridView_(gridView), elementMap_(gridView) {}
- /* \brief Get the finite volume geometry of an element
- * \note The finite volume geometry offers iterators over the sub control volumes
- * and the sub control volume faces of an element.
- */
- const FVElementGeometry& fvGeometry(const Element& element) const
- {
- return fvGeometries_[problem_().elementMapper().index(element)];
- }
-
- //! Get a sub control volume with a global scv index
- const SubControlVolume& subControlVolume(IndexType scvIdx) const
- {
- return *scvs_[scvIdx];
- }
-
- //! Get a sub control volume face with a global scvf index
- const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
- {
- return scvfs_[scvfIdx];
- }
-
//! The total number of sub control volumes
std::size_t numScv() const
{
@@ -105,6 +88,10 @@ public:
Element element(IndexType eIdx) const
{ return elementMap_.element(eIdx); }
+ //! Return the gridView this global object lives on
+ const GridView& gridView() const
+ { return gridView_; }
+
//! update all fvElementGeometries (do this again after grid adaption)
void update(const Problem& problem)
{
@@ -113,7 +100,7 @@ public:
// clear containers (necessary after grid refinement)
scvs_.clear();
scvfs_.clear();
- fvGeometries_.clear();
+ scvfIndicesOfScv_.clear();
elementMap_.clear();
// determine size of containers
@@ -124,8 +111,9 @@ public:
// reserve memory
elementMap_.resize(numScvs);
- scvs_.resize(numScvs);
+ scvs_.reserve(numScvs);
scvfs_.reserve(numScvf);
+ scvfIndicesOfScv_.resize(numScvs);
// Build the scvs and scv faces
IndexType scvfIdx = 0;
@@ -133,7 +121,7 @@ public:
for (const auto& element : elements(gridView_))
{
auto eIdx = problem.elementMapper().index(element);
- scvs_[eIdx] = std::make_shared<SubControlVolume>(element.geometry(), eIdx);
+ scvs_.emplace_back(element.geometry(), eIdx);
// fill the element map with seeds
elementMap_[eIdx] = element.seed();
@@ -169,17 +157,41 @@ public:
}
}
- // Compute the finite volume element geometries
- fvGeometries_.push_back(FVElementGeometry(*this, {eIdx}, scvfsIndexSet));
+ // Save the scvf indices belonging to this scv to build up fv element geometries fast
+ scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
}
}
- // Binding of an element, called by the local jacobian to prepare element assembly
- // For compatibility reasons with the FVGeometry cache disabled
- void bind(const Element& element) {}
- void bindElement(const Element& element) {}
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend FVElementGeometry localView(const CCTpfaGlobalFVGeometry& global)
+ {
+ return FVElementGeometry(global);
+ }
private:
+
+ //! Get a sub control volume with a global scv index
+ const SubControlVolume& scv(IndexType scvIdx) const
+ {
+ return scvs_[scvIdx];
+ }
+
+ //! Get a sub control volume face with a global scvf index
+ const SubControlVolumeFace& scvf(IndexType scvfIdx) const
+ {
+ return scvfs_[scvfIdx];
+ }
+
+ //! Get the sub control volume face indices of an scv by global index
+ const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
+ {
+ return scvfIndicesOfScv_[scvIdx];
+ }
+
const Problem& problem_() const
{ return *problemPtr_; }
@@ -187,15 +199,15 @@ private:
GridView gridView_;
Dumux::ElementMap<GridView> elementMap_;
- std::vector<std::shared_ptr<SubControlVolume>> scvs_;
+ std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
- std::vector<FVElementGeometry> fvGeometries_;
+ std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
IndexType numBoundaryScvf_;
};
// specialization in case the FVElementGeometries are not stored
template<class TypeTag>
-class CCTpfaFVElementGeometryVector<TypeTag, false>
+class CCTpfaGlobalFVGeometry<TypeTag, false>
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -204,32 +216,13 @@ class CCTpfaFVElementGeometryVector<TypeTag, false>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
+ //! The local fvGeometry needs access to the problem
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
public:
//! Constructor
- CCTpfaFVElementGeometryVector(const GridView gridView)
- : gridView_(gridView), elementMap_(gridView), eIdxBound_(-1) {}
-
- /* \brief Get the finite volume geometry of an element
- * \note The finite volume geometry offers iterators over the sub control volumes
- * and the sub control volume faces of an element.
- */
- const FVElementGeometry& fvGeometry(const Element& element) const
- {
- return fvGeometries_[getLocalScvIdx_(problem_().elementMapper().index(element), true)];
- }
-
- //! Get a sub control volume with a global scv index
- const SubControlVolume& subControlVolume(IndexType scvIdx) const
- {
- return localSvs_[getLocalScvIdx_(scvIdx)];
- }
-
- //! Get a sub control volume face with a global scvf index
- const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
- {
- return localScvfs_[getLocalScvFaceIdx_(scvfIdx)];
- }
+ CCTpfaGlobalFVGeometry(const GridView gridView)
+ : gridView_(gridView), elementMap_(gridView) {}
//! The total number of sub control volumes
std::size_t numScv() const
@@ -257,19 +250,22 @@ public:
Element element(IndexType eIdx) const
{ return elementMap_.element(eIdx); }
+ //! Return the gridView this global object lives on
+ const GridView& gridView() const
+ { return gridView_; }
+
//! update all fvElementGeometries (do this again after grid adaption)
void update(const Problem& problem)
{
problemPtr_ = &problem;
elementMap_.clear();
- eIdxBound_ = -1;
// reserve memory or resize the containers
numScvs_ = gridView_.size(0);
numScvf_ = 0;
numBoundaryScvf_ = 0;
elementMap_.resize(numScvs_);
- scvFaceIndices_.resize(numScvs_);
+ scvfIndicesOfScv_.resize(numScvs_);
neighborVolVarIndices_.resize(numScvs_);
// Build the SCV and SCV face
@@ -303,113 +299,33 @@ public:
}
// store the sets of indices in the data container
- scvFaceIndices_[eIdx] = scvfsIndexSet;
+ scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
}
}
- // Binding of an element preparing the geometries of the whole stencil
- // called by the local jacobian to prepare element assembly
- void bind(const Element& element)
- {
- eIdxBound_ = problem_().elementMapper().index(element);
- release_();
- makeElementGeometries_(element);
- for (const auto& intersection : intersections(gridView_, element))
- {
- if (intersection.neighbor())
- makeElementGeometries_(intersection.outside());
- }
- }
-
- // Binding of an element preparing the geometries only inside the element
- void bindElement(const Element& element)
- {
- eIdxBound_ = problem_().elementMapper().index(element);
- release_();
- makeElementGeometries_(element);
- }
-
-
-private:
-
- void release_()
- {
- localSvs_.clear();
- localScvIndices_.clear();
- localScvfs_.clear();
- localScvfIndices_.clear();
- fvGeometries_.clear();
- }
-
- void makeElementGeometries_(const Element& element)
- {
- const auto eIdx = problem_().elementMapper().index(element);
-
- SubControlVolume scv(element.geometry(), eIdx);
- localSvs_.push_back(scv);
- localScvIndices_.push_back(eIdx);
-
- const auto& scvFaceIndices = scvFaceIndices_[eIdx];
- const auto& neighborVolVarIndices = neighborVolVarIndices_[eIdx];
-
- int scvfCounter = 0;
- for (const auto& intersection : intersections(gridView_, element))
- {
- if (intersection.neighbor() || intersection.boundary())
- {
- localScvfs_.push_back(SubControlVolumeFace(intersection.geometry(),
- intersection.geometry().center(),
- intersection.centerUnitOuterNormal(),
- scvFaceIndices[scvfCounter],
- intersection.indexInInside(),
- std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}),
- intersection.boundary()));
-
- localScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
- scvfCounter++;
- }
- }
-
- // Compute the finite volume element geometries
- fvGeometries_.push_back(FVElementGeometry(*this, {eIdx}, scvFaceIndices));
- }
-
- const int getLocalScvIdx_(const int scvIdx, bool isElement = false) const
- {
- auto it = std::find(localScvIndices_.begin(), localScvIndices_.end(), scvIdx);
+ const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
+ { return scvfIndicesOfScv_[scvIdx]; }
- if (it != localScvIndices_.end())
- return std::distance(localScvIndices_.begin(), it);
- else
- {
- if (!isElement)
- DUNE_THROW(Dune::InvalidStateException,
- "Could not find the sub control volume for scvIdx = " << scvIdx <<
- ", make sure to properly bind the FVGeometries to the element before using them.");
- else
- DUNE_THROW(Dune::InvalidStateException,
- "You are trying to get the fvgeometry of element with index " << scvIdx <<
- ", but bound is the element " << eIdxBound_ << ". Please bind element before using the fvgeometry.");
- }
- }
+ const std::vector<IndexType>& neighborVolVarIndices(IndexType scvIdx) const
+ { return neighborVolVarIndices_[scvIdx]; }
- const int getLocalScvFaceIdx_(const int scvfIdx) const
+ /*!
+ * \brief Return a local restriction of this global object
+ * The local object is only functional after calling its bind/bindElement method
+ * This is a free function that will be found by means of ADL
+ */
+ friend FVElementGeometry localView(const CCTpfaGlobalFVGeometry& global)
{
- auto it = std::find(localScvfIndices_.begin(), localScvfIndices_.end(), scvfIdx);
-
- if (it != localScvfIndices_.end())
- return std::distance(localScvfIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException,
- "Could not find the sub control volume face for scvfIdx = " << scvfIdx <<
- ", make sure to properly bind the FVGeometries to the to the element before using them.");
+ return FVElementGeometry(global);
}
+private:
const Problem& problem_() const
{ return *problemPtr_; }
const Problem* problemPtr_;
+
GridView gridView_;
// Information on the global number of geometries
@@ -419,16 +335,8 @@ private:
// vectors that store the global data
Dumux::ElementMap<GridView> elementMap_;
- std::vector<std::vector<IndexType>> scvFaceIndices_;
+ std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
std::vector<std::vector<IndexType>> neighborVolVarIndices_;
-
- // vectors to store the geometries temporarily after binding an element
- IndexType eIdxBound_;
- std::vector<SubControlVolume> localSvs_;
- std::vector<IndexType> localScvIndices_;
- std::vector<SubControlVolumeFace> localScvfs_;
- std::vector<IndexType> localScvfIndices_;
- std::vector<FVElementGeometry> fvGeometries_;
};
} // end namespace
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh
index aaba6ee555..40d78cfaee 100644
--- a/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh
@@ -43,9 +43,14 @@ private:
public:
// the contructor in the cc case
- CCTpfaSubControlVolume(Geometry geometry,
- IndexType elementIndex)
- : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)) {}
+ CCTpfaSubControlVolume(Geometry&& geometry,
+ IndexType elementIndex)
+ : SubControlVolumeBase<G, I>(std::move(geometry), elementIndex) {}
+
+ // the contructor in the cc case
+ CCTpfaSubControlVolume(const Geometry& geometry,
+ IndexType elementIndex)
+ : SubControlVolumeBase<G, I>(geometry, elementIndex) {}
IndexType indexInElement() const
{
diff --git a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
index 9a3d763c5c..6c43bc5015 100644
--- a/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
+++ b/dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh
@@ -38,6 +38,7 @@ namespace Dumux
template<class Geometry, typename IndexType>
class CCTpfaSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, IndexType>
{
+ using ParentType = SubControlVolumeFaceBase<Geometry, IndexType>;
using Scalar = typename Geometry::ctype;
static const int dim = Geometry::mydimension;
static const int dimworld = Geometry::coorddimension;
@@ -46,14 +47,27 @@ class CCTpfaSubControlVolumeFace : public SubControlVolumeFaceBase<Geometry, Ind
using LocalPosition = Dune::FieldVector<Scalar, dim>;
public:
+ //! Constructor with all arguments
CCTpfaSubControlVolumeFace(const Geometry& geometry,
const GlobalPosition& ipGlobal,
const GlobalPosition& unitOuterNormal,
IndexType scvfIndex,
- IndexType indexInElement,
const std::vector<IndexType>& scvIndices,
bool boundary = false)
- : SubControlVolumeFaceBase<Geometry, IndexType>(geometry, ipGlobal, unitOuterNormal, scvfIndex, indexInElement, scvIndices, boundary)
+ : ParentType(geometry, ipGlobal, unitOuterNormal, scvfIndex, scvIndices, boundary)
+ {}
+
+ //! Constructor with intersection
+ template <class Intersection>
+ CCTpfaSubControlVolumeFace(const Intersection& is,
+ IndexType scvfIndex,
+ const std::vector<IndexType>& scvIndices)
+ : ParentType(is.geometry(),
+ is.geometry().center(),
+ is.centerUnitOuterNormal(),
+ scvfIndex,
+ scvIndices,
+ is.boundary())
{}
};
diff --git a/dumux/discretization/cellcentered/volumevariablesvector.hh b/dumux/discretization/cellcentered/volumevariablesvector.hh
index 39b6235652..1fcc38cf14 100644
--- a/dumux/discretization/cellcentered/volumevariablesvector.hh
+++ b/dumux/discretization/cellcentered/volumevariablesvector.hh
@@ -78,7 +78,7 @@ public:
{
problem.model().fvGeometries_().bindElement(element);
const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
(*this)[scv].update(sol[scv.dofIndex()],
problem,
@@ -87,7 +87,7 @@ public:
}
// handle the boundary volume variables
- for (const auto& scvFace : scvfs(fvGeometry))
+ for (auto&& scvFace : scvfs(fvGeometry))
{
// if we are not on a boundary, skip the rest
if (!scvFace.boundary())
@@ -161,6 +161,7 @@ class CCVolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
static const int dim = GridView::dimension;
@@ -187,7 +188,8 @@ public:
// Binding of an element, prepares the volume variables within the element stencil
// called by the local jacobian to prepare element assembly
- void bind(const Element& element)
+ void bind(const Element& element,
+ const FVElementGeometry& fvGeometry)
{
eIdxBound_ = problem_().elementMapper().index(element);
@@ -202,7 +204,7 @@ public:
int localIdx = 0;
// update the volume variables of the element at hand
- const auto& scvI = problem_().model().fvGeometries().subControlVolume(eIdxBound_);
+ auto&& scvI = fvGeometry.scv(eIdxBound_);
const auto& solI = problem_().model().curSol()[eIdxBound_];
volumeVariables_[localIdx].update(solI, problem_(), element, scvI);
volVarIndices_[localIdx] = scvI.index();
@@ -211,8 +213,8 @@ public:
// Update the volume variables of the neighboring elements
for (auto globalJ : neighborStencil)
{
- const auto& elementJ = problem_().model().fvGeometries().element(globalJ);
- const auto& scvJ = problem_().model().fvGeometries().subControlVolume(globalJ);
+ const auto& elementJ = fvGeometry.globalFvGeometry().element(globalJ);
+ auto&& scvJ = fvGeometry.scv(globalJ);
const auto& solJ = problem_().model().curSol()[globalJ];
volumeVariables_[localIdx].update(solJ, problem_(), elementJ, scvJ);
volVarIndices_[localIdx] = scvJ.index();
@@ -220,8 +222,7 @@ public:
}
// Update boundary volume variables
- const auto& fvGeometry = problem_().model().fvGeometries(element);
- for (const auto& scvFace : scvfs(fvGeometry))
+ for (auto&& scvFace : scvfs(fvGeometry))
{
// if we are not on a boundary, skip the rest
if (!scvFace.boundary())
@@ -243,14 +244,15 @@ public:
// Binding of an element, prepares only the volume variables of the element
// specialization for cc models
- void bindElement(const Element& element)
+ void bindElement(const Element& element,
+ const FVElementGeometry& fvGeometry)
{
eIdxBound_ = problem_().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element at hand
- const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdxBound_);
+ auto&& scv = fvGeometry.scv(eIdxBound_);
const auto& sol = problem_().model().curSol()[eIdxBound_];
volumeVariables_[0].update(sol, problem_(), element, scv);
volVarIndices_[0] = scv.index();
@@ -283,12 +285,8 @@ private:
const int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
-
- if (it != volVarIndices_.end())
- return std::distance(volVarIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not find the current volume variables for volVarIdx = " << volVarIdx <<
- ", make sure to properly bind the volume variables to the element before using them");
+ assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
+ return std::distance(volVarIndices_.begin(), it);
}
Problem& problem_() const
@@ -312,6 +310,7 @@ class CCVolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using IndexType = typename GridView::IndexSet::IndexType;
static const int dim = GridView::dimension;
@@ -337,14 +336,15 @@ public:
}
// Binding of an element, prepares the volume variables of only the element
- void bindElement(const Element& element)
+ void bindElement(const Element& element,
+ const FVElementGeometry& fvGeometry)
{
eIdxBound_ = problem_().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element at hand
- const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdxBound_);
+ auto&& scv = fvGeometry.scv(eIdxBound_);
const auto& sol = problem_().model().prevSol()[eIdxBound_];
volumeVariables_[0].update(sol, problem_(), element, scv);
volVarIndices_[0] = scv.index();
@@ -377,12 +377,8 @@ private:
const int getLocalIdx_(const int volVarIdx) const
{
auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
-
- if (it != volVarIndices_.end())
- return std::distance(volVarIndices_.begin(), it);
- else
- DUNE_THROW(Dune::InvalidStateException, "Could not find the previous volume variables for volVarIdx = " << volVarIdx <<
- ", make sure to properly bind the volume variables to the element before using them");
+ assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
+ return std::distance(volVarIndices_.begin(), it);
}
Problem& problem_() const
diff --git a/dumux/discretization/fluxvariablesbase.hh b/dumux/discretization/fluxvariablesbase.hh
index de8061254d..52eb6dfd72 100644
--- a/dumux/discretization/fluxvariablesbase.hh
+++ b/dumux/discretization/fluxvariablesbase.hh
@@ -42,6 +42,7 @@ class FluxVariablesBase
using IndexType = typename GridView::IndexSet::IndexType;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Stencil = std::vector<IndexType>;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
enum{ enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache) };
@@ -53,15 +54,17 @@ public:
void init(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace &scvFace)
{
problemPtr_ = &problem;
elementPtr_ = &element;
scvFacePtr_ = &scvFace;
+ fvGeometryPtr_ = &fvGeometry;
// update the stencil if needed
if (!enableFluxVarsCache)
- stencil_ = asImp_().computeStencil(problem, element, scvFace);
+ stencil_ = asImp_().computeStencil(problem, element, fvGeometry, scvFace);
}
// when caching is enabled, get the stencil from the cache class
@@ -83,6 +86,9 @@ public:
const SubControlVolumeFace& scvFace() const
{ return *scvFacePtr_; }
+ const FVElementGeometry& fvGeometry() const
+ { return *fvGeometryPtr_; }
+
Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
{ DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); }
@@ -101,7 +107,8 @@ private:
}
const Problem* problemPtr_; //! Pointer to the problem
- const Element* elementPtr_; //! Pointer to the element at hand
+ const Element* elementPtr_; //! Pointer to the element at hand
+ const FVElementGeometry* fvGeometryPtr_;
const SubControlVolumeFace* scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created
Stencil stencil_; //! The flux stencil
};
diff --git a/dumux/discretization/fvelementgeometry.hh b/dumux/discretization/fvelementgeometry.hh
deleted file mode 100644
index 7640c0be3f..0000000000
--- a/dumux/discretization/fvelementgeometry.hh
+++ /dev/null
@@ -1,207 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief Class providing iterators over sub control volumes and sub control
- * volume faces of an element.
- */
-#ifndef DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH
-#define DUMUX_DISCRETIZATION_FV_ELEMENTGEOMETRY_HH
-
-#include <dune/common/iteratorrange.hh>
-#include <dune/common/iteratorfacades.hh>
-
-namespace Dumux
-{
-namespace Properties
-{
-NEW_PROP_TAG(SubControlVolume);
-NEW_PROP_TAG(SubControlVolumeFace);
-NEW_PROP_TAG(FVElementGeometry);
-NEW_PROP_TAG(FVElementGeometryVector);
-}
-
-/*!
- * \ingroup Discretization
- * \brief An iterator over sub control volumes
- */
-template<class SubControlVolume, class Vector, class FVElementGeometryVector>
-class ScvIterator : public Dune::ForwardIteratorFacade<ScvIterator<SubControlVolume,
- Vector,
- FVElementGeometryVector>,
- const SubControlVolume>
-{
- using ThisType = ScvIterator<SubControlVolume, Vector, FVElementGeometryVector>;
- using Iterator = typename Vector::const_iterator;
-public:
- ScvIterator(const Iterator& it, const FVElementGeometryVector& fvGeometryVector)
- : it_(it), fvGeometryVector_(&fvGeometryVector) {}
-
- //! default constructor
- ScvIterator() : it_(Iterator()), fvGeometryVector_(nullptr) {}
-
- const SubControlVolume& dereference() const
- {
- return fvGeometryVector_->subControlVolume(*it_);
- }
-
- bool equals(const ThisType& other) const
- {
- return it_ == other.it_;
- }
-
- void increment()
- {
- ++it_;
- }
-
-private:
- Iterator it_;
- const FVElementGeometryVector* fvGeometryVector_;
-};
-
-/*!
- * \ingroup ImplcititModel
- * \brief An iterator over sub control volume faces
- */
-template<class SubControlVolumeFace, class Vector, class FVElementGeometryVector>
-class ScvfIterator : public Dune::ForwardIteratorFacade<ScvfIterator<SubControlVolumeFace,
- Vector,
- FVElementGeometryVector>,
- const SubControlVolumeFace>
-{
- using ThisType = ScvfIterator<SubControlVolumeFace, Vector, FVElementGeometryVector>;
- using Iterator = typename Vector::const_iterator;
-public:
- ScvfIterator(const Iterator& it, const FVElementGeometryVector& fvGeometryVector)
- : it_(it), fvGeometryVector_(&fvGeometryVector) {}
-
- //! default constructor
- ScvfIterator() : it_(Iterator()), fvGeometryVector_(nullptr) {}
-
- const SubControlVolumeFace& dereference() const
- {
- return fvGeometryVector_->subControlVolumeFace(*it_);
- }
-
- bool equals(const ThisType& other) const
- {
- return it_ == other.it_;
- }
-
- void increment()
- {
- it_++;
- }
-
-private:
- Iterator it_;
- const FVElementGeometryVector* fvGeometryVector_;
-};
-
-/*!
- * \ingroup ImplcititModel
- * \brief Provide iterators over sub control volumes and sub control
- * volume faces of an element..
- */
-template<class TypeTag>
-class FVElementGeometry
-{
- using ThisType = Dumux::FVElementGeometry<TypeTag>;
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVElementGeometryVector = typename GET_PROP_TYPE(TypeTag, FVElementGeometryVector);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, FVElementGeometryVector>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, FVElementGeometryVector>;
-
-public:
- // This class in not default-constructible
- FVElementGeometry() = delete;
-
- // Constructor with vectors
- FVElementGeometry(const FVElementGeometryVector& localFvGeometry,
- const std::vector<IndexType>& scvIndices,
- const std::vector<IndexType>& scvfIndices)
- : localFvGeometry_(localFvGeometry),
- scvIndices_(scvIndices),
- scvfIndices_(scvfIndices)
- {}
-
- //! iterator range for sub control volumes
- //! This is a free function found by means of ADL
- //! To iterate over all sub control volumes of this FVElementGeometry use
- //! for (const auto& scv : scvs(fvGeometry))
- friend inline Dune::IteratorRange<ScvIterator>
- scvs(const FVElementGeometry& g)
- {
- return Dune::IteratorRange<ScvIterator>(ScvIterator(g.scvIndices().begin(), g.localFvGeometry()),
- ScvIterator(g.scvIndices().end(), g.localFvGeometry()));
- }
-
- //! iterator range for sub control volume faces
- //! This is a free function found by means of ADL
- //! To iterate over all sub control volume faces of this FVElementGeometry use
- //! for (const auto& scvf : scvfs(fvGeometry))
- friend inline Dune::IteratorRange<ScvfIterator>
- scvfs(const FVElementGeometry& g)
- {
- return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndices().begin(), g.localFvGeometry()),
- ScvfIterator(g.scvfIndices().end(), g.localFvGeometry()));
- }
-
- //! number of sub control volumes in this fv element geometry
- std::size_t numScv() const
- {
- return scvIndices_.size();
- }
-
- //! number of sub control volumes in this fv element geometry
- std::size_t numScvf() const
- {
- return scvfIndices_.size();
- }
-
-private:
-
- //! The sub control volume indices
- //! Depending on the type of discretization and caching
- //! these may be local or global indices
- const std::vector<IndexType>& scvIndices() const
- { return scvIndices_; }
-
- //! The sub control volume face indices
- //! Depending on the type of discretization and caching
- //! these may be local or global indices
- const std::vector<IndexType>& scvfIndices() const
- { return scvfIndices_; }
-
- //! The LocalFvGeometry object this FvElementGeometry belongs to
- const FVElementGeometryVector& localFvGeometry() const
- { return localFvGeometry_; }
-
- const FVElementGeometryVector& localFvGeometry_;
- std::vector<IndexType> scvIndices_;
- std::vector<IndexType> scvfIndices_;
-};
-
-}
-
-#endif
diff --git a/dumux/discretization/scvandscvfiterators.hh b/dumux/discretization/scvandscvfiterators.hh
new file mode 100644
index 0000000000..81c01a1ab8
--- /dev/null
+++ b/dumux/discretization/scvandscvfiterators.hh
@@ -0,0 +1,112 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Class providing iterators over sub control volumes and sub control
+ * volume faces of an element.
+ */
+#ifndef DUMUX_SCV_AND_SCVF_ITERATORS_HH
+#define DUMUX_SCV_AND_SCVF_ITERATORS_HH
+
+#include <dune/common/iteratorfacades.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Discretization
+ * \brief An iterator over sub control volumes
+ */
+template<class SubControlVolume, class Vector, class FVElementGeometry>
+class ScvIterator : public Dune::ForwardIteratorFacade<ScvIterator<SubControlVolume,
+ Vector,
+ FVElementGeometry>,
+ const SubControlVolume>
+{
+ using ThisType = ScvIterator<SubControlVolume, Vector, FVElementGeometry>;
+ using Iterator = typename Vector::const_iterator;
+public:
+ ScvIterator(const Iterator& it, const FVElementGeometry& fvGeometry)
+ : it_(it), fvGeometryPtr_(&fvGeometry) {}
+
+ //! default constructor
+ ScvIterator() : it_(Iterator()), fvGeometryPtr_(nullptr) {}
+
+ const SubControlVolume& dereference() const
+ {
+ return fvGeometryPtr_->scv(*it_);
+ }
+
+ bool equals(const ThisType& other) const
+ {
+ return it_ == other.it_;
+ }
+
+ void increment()
+ {
+ ++it_;
+ }
+
+private:
+ Iterator it_;
+ const FVElementGeometry* fvGeometryPtr_;
+};
+
+/*!
+ * \ingroup ImplcititModel
+ * \brief An iterator over sub control volume faces
+ */
+template<class SubControlVolumeFace, class Vector, class FVElementGeometry>
+class ScvfIterator : public Dune::ForwardIteratorFacade<ScvfIterator<SubControlVolumeFace,
+ Vector,
+ FVElementGeometry>,
+ const SubControlVolumeFace>
+{
+ using ThisType = ScvfIterator<SubControlVolumeFace, Vector, FVElementGeometry>;
+ using Iterator = typename Vector::const_iterator;
+public:
+ ScvfIterator(const Iterator& it, const FVElementGeometry& fvGeometry)
+ : it_(it), fvGeometryPtr_(&fvGeometry) {}
+
+ //! default constructor
+ ScvfIterator() : it_(Iterator()), fvGeometryPtr_(nullptr) {}
+
+ const SubControlVolumeFace& dereference() const
+ {
+ return fvGeometryPtr_->scvf(*it_);
+ }
+
+ bool equals(const ThisType& other) const
+ {
+ return it_ == other.it_;
+ }
+
+ void increment()
+ {
+ it_++;
+ }
+
+private:
+ Iterator it_;
+ const FVElementGeometry* fvGeometryPtr_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/subcontrolvolumebase.hh b/dumux/discretization/subcontrolvolumebase.hh
index 56f141514a..15722ac3e3 100644
--- a/dumux/discretization/subcontrolvolumebase.hh
+++ b/dumux/discretization/subcontrolvolumebase.hh
@@ -32,9 +32,12 @@ namespace Dumux
* \brief Base class for a sub control volume, i.e a part of the control
* volume we are making the balance for.
*/
-template<class Geometry, typename IndexType>
+template<class G, typename I>
class SubControlVolumeBase
{
+ using IndexType = I;
+ using Geometry = typename std::decay<G>::type;
+
using Scalar = typename Geometry::ctype;
enum { dimworld = Geometry::coorddimension };
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
@@ -43,24 +46,54 @@ public:
/*\brief The constructor
* \param geometry The geometry of the sub control volume
* \param elementIndex The index of the element the scv is part of
- * \param indexInElement The local index of the scv in the element
*/
- SubControlVolumeBase(Geometry geometry,
+ SubControlVolumeBase(Geometry&& geometry,
IndexType elementIndex)
: geometry_(std::move(geometry)),
elementIndex_(std::move(elementIndex))
{}
+ /*\brief The constructor
+ * \param geometry The geometry of the sub control volume
+ * \param elementIndex The index of the element the scv is part of
+ */
+ SubControlVolumeBase(const Geometry& geometry,
+ IndexType elementIndex)
+ : geometry_(geometry),
+ elementIndex_(elementIndex)
+ {}
+
+ //! The copy constrcutor
+ SubControlVolumeBase(const SubControlVolumeBase& other) = delete;
+
+ //! The move constrcutor
+ SubControlVolumeBase(SubControlVolumeBase&& other) = default;
+
+ //! The copy assignment operator
+ SubControlVolumeBase& operator=(const SubControlVolumeBase& other) = delete;
+
+ //! The move assignment operator
+ SubControlVolumeBase& operator=(SubControlVolumeBase&& other)
+ {
+ // We want to use the default copy/move assignment operators.
+ // But since geometry is not assignable :( we
+ // have to dirty reconstruct it
+ geometry_.~Geometry();
+ new(&geometry_) Geometry(std::move(other.geometry_));
+ elementIndex_ = std::move(other.elementIndex_);
+ return *this;
+ }
+
//! The center of the sub control volume
GlobalPosition center() const
{
- return geometry_.center();
+ return geometry().center();
}
//! The volume of the sub control volume
Scalar volume() const
{
- return geometry_.volume();
+ return geometry().volume();
}
//! The geometry of the sub control volume
@@ -81,104 +114,6 @@ private:
IndexType elementIndex_;
};
-// This class is set as the property
-template<class G, typename I, bool isBox>
-class SubControlVolume {};
-
-// Specialization for cc models
-template<class G, typename I>
-class SubControlVolume<G, I, false> : public SubControlVolumeBase<G, I>
-{
-public:
- // exported types
- using Geometry = G;
- using IndexType = I;
-
-private:
- using Scalar = typename Geometry::ctype;
- enum { dimworld = Geometry::coorddimension };
- using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
-
-public:
- // the contructor in the cc case
- SubControlVolume(Geometry geometry,
- IndexType elementIndex)
- : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)) {}
-
- IndexType indexInElement() const
- {
- return IndexType(0);
- }
-
- //! The global index of this scv
- IndexType index() const
- {
- return this->elementIndex();
- }
-
- IndexType dofIndex() const
- {
- return this->elementIndex();
- }
-
- GlobalPosition dofPosition() const
- {
- return this->geometry().center();
- }
-};
-
-// Specialization for box models
-template<class G, typename I>
-class SubControlVolume<G, I, true> : public SubControlVolumeBase<G, I>
-{
-public:
- // exported types
- using Geometry = G;
- using IndexType = I;
-
-private:
- using Scalar = typename Geometry::ctype;
- enum { dimworld = Geometry::coorddimension };
- using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
-
-public:
- // the contructor in the box case
- SubControlVolume(Geometry geometry,
- IndexType scvIdx,
- IndexType elementIndex,
- IndexType indexInElement,
- IndexType dofIndex)
- : SubControlVolumeBase<G, I>(std::move(geometry), std::move(elementIndex)),
- scvIdx_(scvIdx), indexInElement_(std::move(indexInElement)),
- dofIndex_(std::move(dofIndex)) {}
-
- IndexType indexInElement() const
- {
- return indexInElement_;
- }
-
- //! The global index of this scv
- IndexType index() const
- {
- return scvIdx_;
- }
-
- IndexType dofIndex() const
- {
- return dofIndex_;
- }
-
- GlobalPosition dofPosition() const
- {
- return this->geometry().corner(indexInElement_);
- }
-
-private:
- IndexType scvIdx_;
- IndexType indexInElement_;
- IndexType dofIndex_;
-};
-
} // end namespace
#endif
diff --git a/dumux/discretization/subcontrolvolumefacebase.hh b/dumux/discretization/subcontrolvolumefacebase.hh
index ab732275d7..01e962c315 100644
--- a/dumux/discretization/subcontrolvolumefacebase.hh
+++ b/dumux/discretization/subcontrolvolumefacebase.hh
@@ -49,7 +49,6 @@ public:
const GlobalPosition& ipGlobal,
const GlobalPosition& unitOuterNormal,
IndexType scvfIndex,
- IndexType indexInElement,
const std::vector<IndexType>& scvIndices,
bool boundary = false)
: geometry_(geometry),
@@ -57,7 +56,20 @@ public:
ipLocal_(geometry.local(ipGlobal)),
unitOuterNormal_(unitOuterNormal),
scvfIndex_(scvfIndex),
- indexInElement_(indexInElement),
+ scvIndices_(scvIndices),
+ boundary_(boundary) {}
+
+ SubControlVolumeFaceBase(Geometry&& geometry,
+ const GlobalPosition& ipGlobal,
+ const GlobalPosition& unitOuterNormal,
+ IndexType scvfIndex,
+ const std::vector<IndexType>& scvIndices,
+ bool boundary = false)
+ : geometry_(std::move(geometry)),
+ ipGlobal_(ipGlobal),
+ ipLocal_(geometry_.local(ipGlobal)),
+ unitOuterNormal_(unitOuterNormal),
+ scvfIndex_(scvfIndex),
scvIndices_(scvIndices),
boundary_(boundary) {}
@@ -123,19 +135,12 @@ public:
return scvfIndex_;
}
- //! The global index of this sub control volume face
- IndexType indexInElement() const
- {
- return indexInElement_;
- }
-
private:
Geometry geometry_;
GlobalPosition ipGlobal_;
LocalPosition ipLocal_;
GlobalPosition unitOuterNormal_;
IndexType scvfIndex_;
- IndexType indexInElement_;
std::vector<IndexType> scvIndices_;
bool boundary_;
};
diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh
index 809b59ffc4..b2c5c13b1b 100644
--- a/dumux/geomechanics/elastic/model.hh
+++ b/dumux/geomechanics/elastic/model.hh
@@ -166,6 +166,7 @@ public:
sigmay[eIdx] += stress[1];
if (dim == 3)
sigmaz[eIdx] += stress[2];
+ }
}
writer.attachDofData(ux, "ux", isBox);
diff --git a/dumux/implicit/assembler.hh b/dumux/implicit/assembler.hh
index 09e16ae56f..e824559e46 100644
--- a/dumux/implicit/assembler.hh
+++ b/dumux/implicit/assembler.hh
@@ -24,6 +24,7 @@
#define DUMUX_IMPLICIT_ASSEMBLER_HH
#include "properties.hh"
+#include <dune/istl/io.hh>
namespace Dumux {
@@ -115,6 +116,8 @@ public:
DUNE_THROW(NumericalProblem,
"A process did not succeed in linearizing the system");
}
+ // printmatrix(std::cout, matrix(), "", "");
+
}
/*!
diff --git a/dumux/implicit/box/elementboundarytypes.hh b/dumux/implicit/box/elementboundarytypes.hh
index 69ec63c1c0..1a8e24d072 100644
--- a/dumux/implicit/box/elementboundarytypes.hh
+++ b/dumux/implicit/box/elementboundarytypes.hh
@@ -89,9 +89,9 @@ public:
hasNeumann_ = false;
hasOutflow_ = false;
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
- int scvIdxLocal = scv.indexInElement();
+ int scvIdxLocal = scv.index();
(*this)[scvIdxLocal].reset();
if (problem.model().onBoundary(scv))
diff --git a/dumux/implicit/box/localjacobian.hh b/dumux/implicit/box/localjacobian.hh
index bd1b8af5bc..38ef7f2044 100644
--- a/dumux/implicit/box/localjacobian.hh
+++ b/dumux/implicit/box/localjacobian.hh
@@ -68,30 +68,30 @@ namespace Dumux
template<class TypeTag>
class BoxLocalJacobian : public ImplicitLocalJacobian<TypeTag>
{
- typedef ImplicitLocalJacobian<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
+ using ParentType = ImplicitLocalJacobian<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, LocalJacobian);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Model = typename GET_PROP_TYPE(TypeTag, Model);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using JacobianAssembler = typename GET_PROP_TYPE(TypeTag, JacobianAssembler);
enum {
numEq = GET_PROP_VALUE(TypeTag, NumEq),
dim = GridView::dimension,
};
- typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
- typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
BoxLocalJacobian()
@@ -105,41 +105,41 @@ public:
*
* \param element The DUNE Codim<0> entity which we look at.
*/
- void assemble(const Element& element, JacobianMatrix& matrix, SolutionVector& residual)
+ void assemble(const Element& element,
+ JacobianMatrix& matrix,
+ SolutionVector& residual)
{
// prepare the volvars/fvGeometries for the case when caching is disabled
- this->model_().fvGeometries_().bind(element);
- this->model_().curVolVars_().bind(element);
- this->model_().prevVolVars_().bindElement(element);
- this->model_().fluxVariablesCache_().bind(element);
-
- // the finite volume element geometry
- const auto& fvGeometry = this->model_().fvGeometries(element);
+ auto fvGeometry = localView(this->problem_().model().globalFvGeometry());
+ fvGeometry.bind(element);
+ this->model_().curVolVars_().bind(element, fvGeometry);
+ this->model_().prevVolVars_().bindElement(element, fvGeometry);
+ this->model_().fluxVariablesCache_().bind(element, fvGeometry);
// the element boundary types
- bcTypes_.update(this->problem_(), element, fvGeometry);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeometry);
// calculate the actual element residual
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes, fvGeometry);
this->residual_ = this->localResidual().residual();
this->model_().updatePVWeights(fvGeometry);
// calculation of the derivatives
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
// dof index and corresponding actual pri vars
const auto dofIdx = scv.dofIndex();
- const auto localScvIdx = scv.indexInElement();
VolumeVariables origVolVars(this->model_().curVolVars(scv));
// add precalculated residual for this scv into the global container
- residual[dofIdx] += this->residual_[localScvIdx];
+ residual[dofIdx] += this->residual_[scv.index()];
// calculate derivatives w.r.t to the privars at the dof at hand
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
{
- evalPartialDerivative_(matrix, element, scv, pvIdx);
+ evalPartialDerivative_(matrix, element, fvGeometry, scv, bcTypes, pvIdx);
// restore the original state of the scv's volume variables
this->model_().curVolVars_(scv) = origVolVars;
@@ -195,10 +195,12 @@ protected:
*/
void evalPartialDerivative_(JacobianMatrix& matrix,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
+ const ElementBoundaryTypes& bcTypes,
const int pvIdx)
{
- auto dofIdx = scv.dofIndex();
+ const auto dofIdx = scv.dofIndex();
ElementSolutionVector partialDeriv(element.subEntities(dim));
PrimaryVariables priVars(this->model_().curSol()[dofIdx]);
@@ -219,7 +221,7 @@ protected:
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
// calculate the deflected residual
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes, fvGeometry);
// store the residual
partialDeriv = this->localResidual().residual();
@@ -245,7 +247,7 @@ protected:
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
// calculate the deflected residual
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes, fvGeometry);
// subtract the residual from the derivative storage
partialDeriv -= this->localResidual().residual();
@@ -263,17 +265,14 @@ protected:
partialDeriv /= delta;
// update the global stiffness matrix with the current partial derivatives
- const auto& fvGeometry = this->model_().fvGeometries(element);
- for (const auto& scvJ : scvs(fvGeometry))
- this->updateGlobalJacobian_(matrix, scvJ.dofIndex(), dofIdx, pvIdx, partialDeriv[scvJ.indexInElement()]);
+ for (auto&& scvJ : scvs(fvGeometry))
+ this->updateGlobalJacobian_(matrix, scvJ.dofIndex(), dofIdx, pvIdx, partialDeriv[scvJ.index()]);
}
private:
-
int numericDifferenceMethod_;
-
- ElementBoundaryTypes bcTypes_;
};
+
}
#endif
diff --git a/dumux/implicit/box/localresidual.hh b/dumux/implicit/box/localresidual.hh
index 0f20ee7e1c..d8ff006582 100644
--- a/dumux/implicit/box/localresidual.hh
+++ b/dumux/implicit/box/localresidual.hh
@@ -44,45 +44,78 @@ namespace Dumux
template<class TypeTag>
class BoxLocalResidual : public ImplicitLocalResidual<TypeTag>
{
- typedef ImplicitLocalResidual<TypeTag> ParentType;
+ using ParentType = ImplicitLocalResidual<TypeTag>;
friend class ImplicitLocalResidual<TypeTag>;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
enum {
numEq = GET_PROP_VALUE(TypeTag, NumEq),
dim = GridView::dimension
};
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+public:
// copying the local residual class is not a good idea
- BoxLocalResidual(const BoxLocalResidual &);
+ BoxLocalResidual(const BoxLocalResidual &) = delete;
-public:
- BoxLocalResidual() : ParentType()
- { }
+ BoxLocalResidual() : ParentType() {}
protected:
- void evalFluxes_()
+ void evalFluxes_(const Element& element,
+ const ElementBoundaryTypes& bcTypes,
+ const FVElementGeometry& fvGeometry)
{
// calculate the mass flux over the scv faces and subtract
- const auto& fvGeometry = this->fvGeometry_();
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
if (!scvf.boundary())
{
- auto flux = this->asImp_().computeFlux(scvf);
- const auto& insideScv = this->model_().fvGeometries().subControlVolume(scvf.insideScvIdx());
- const auto& outsideScv = this->model_().fvGeometries().subControlVolume(scvf.outsideScvIdx());
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+
+ auto flux = this->asImp_().computeFlux(element, fvGeometry, scvf, bcTypes[insideScv.index()]);
- this->residual_[insideScv.indexInElement()] += flux;
- this->residual_[outsideScv.indexInElement()] -= flux;
+ this->residual_[insideScv.index()] += flux;
+ this->residual_[outsideScv.index()] -= flux;
+ }
+ }
+ }
+
+ void evalBoundary_(const Element &element,
+ const ElementBoundaryTypes& bcTypes,
+ const FVElementGeometry& fvGeometry)
+ {
+ if (bcTypes.hasNeumann() || bcTypes.hasOutflow())
+ {
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (scvf.boundary())
+ {
+ auto&& scv = fvGeometry.scv(scvf.insideScvIdx());
+ this->asImp_().evalBoundaryFluxes_(element, scvf, scv, bcTypes[scv.index()]);
+ }
+ }
+ }
+
+ // additionally treat mixed D/N conditions in a strong sense
+ if (bcTypes.hasDirichlet())
+ {
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ auto scvBcTypes = bcTypes[scv.index()];
+ if (!scvBcTypes.hasDirichlet())
+ continue;
+
+ this->asImp_().evalDirichlet_(element, scv, scvBcTypes);
}
}
}
@@ -91,9 +124,11 @@ protected:
* \brief Set the values of the Dirichlet boundary control volumes
* of the current element.
*/
- void evalDirichlet_(const SubControlVolume& scv, const BoundaryTypes& bcTypes)
+ void evalDirichlet_(const Element& element,
+ const SubControlVolume &scv,
+ const BoundaryTypes& bcTypes)
{
- PrimaryVariables dirichletValues = this->asImp_().problem_().dirichlet(this->element_(), scv);
+ PrimaryVariables dirichletValues = this->problem().dirichlet(element, scv);
// set the dirichlet conditions
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@@ -105,9 +140,9 @@ protected:
Valgrind::CheckDefined(dirichletValues[pvIdx]);
// get the primary variables
- const auto& priVars = this->model_().curVolVars(scv).priVars();
+ const auto& priVars = this->problem().model().curVolVars(scv).priVars();
- this->residual_[scv.indexInElement()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+ this->residual_[scv.index()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
}
}
}
@@ -115,12 +150,15 @@ protected:
/*!
* \brief Add all fluxes resulting from Neumann and outflow boundary conditions to the local residual.
*/
- void evalBoundaryFluxes_(const SubControlVolumeFace &scvf, const SubControlVolume& insideScv, const BoundaryTypes& bcTypes)
+ void evalBoundaryFluxes_(const Element& element,
+ const SubControlVolumeFace &scvf,
+ const SubControlVolume& insideScv,
+ const BoundaryTypes& bcTypes)
{
// evaluate the Neumann conditions at the boundary face
if (bcTypes.hasNeumann())
- this->residual_[insideScv.indexInElement()] += this->asImp_().evalNeumannSegment_(scvf, insideScv, bcTypes);
+ this->residual_[insideScv.index()] += this->asImp_().evalNeumannSegment_(element, scvf, insideScv, bcTypes);
// TODO: evaluate the outflow conditions at the boundary face
//if (bcTypes.hasOutflow())
@@ -130,17 +168,18 @@ protected:
/*!
* \brief Add Neumann boundary conditions for a single scv face
*/
- PrimaryVariables evalNeumannSegment_(const SubControlVolumeFace &scvf,
+ PrimaryVariables evalNeumannSegment_(const Element& element,
+ const SubControlVolumeFace &scvf,
const SubControlVolume& insideScv,
const BoundaryTypes &bcTypes)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables flux(0);
- auto neumannFluxes = this->problem_().neumann(this->element_(), scvf);
+ auto neumannFluxes = this->problem().neumann(element, scvf);
// multiply neumann fluxes with the area and the extrusion factor
- neumannFluxes *= scvf.area()*this->problem_().model().curVolVars(insideScv).extrusionFactor();
+ neumannFluxes *= scvf.area()*this->problem().model().curVolVars(insideScv).extrusionFactor();
// add fluxes to the temporary vector
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@@ -181,35 +220,6 @@ protected:
// }
// }
// }
-
- void evalBoundary_()
- {
- const auto& fvGeometry = this->fvGeometry_();
- if (this->bcTypes_().hasNeumann() || this->bcTypes_().hasOutflow())
- {
- for (const auto& scvf : scvfs(fvGeometry))
- {
- if (scvf.boundary())
- {
- const auto& scv = this->problem_().model().fvGeometries().subControlVolume(scvf.insideScvIdx());
- this->asImp_().evalBoundaryFluxes_(scvf, scv, this->bcTypes_(scv.indexInElement()));
- }
- }
- }
-
- // additionally treat mixed D/N conditions in a strong sense
- if (this->bcTypes_().hasDirichlet())
- {
- for (const auto& scv : scvs(fvGeometry))
- {
- BoundaryTypes bcTypes = this->bcTypes_(scv.indexInElement());
- if (!bcTypes.hasDirichlet())
- continue;
-
- this->asImp_().evalDirichlet_(scv, bcTypes);
- }
- }
- }
};
}
diff --git a/dumux/implicit/box/propertydefaults.hh b/dumux/implicit/box/propertydefaults.hh
index 3681ef036c..71d919fcea 100644
--- a/dumux/implicit/box/propertydefaults.hh
+++ b/dumux/implicit/box/propertydefaults.hh
@@ -27,12 +27,12 @@
#define DUMUX_BOX_PROPERTY_DEFAULTS_HH
#include <dumux/implicit/propertydefaults.hh>
-#include <dumux/discretization/fvelementgeometry.hh>
#include <dumux/discretization/box/subcontrolvolume.hh>
#include <dumux/discretization/box/subcontrolvolumeface.hh>
#include <dumux/discretization/box/fluxvariablescachevector.hh>
#include <dumux/discretization/box/volumevariablesvector.hh>
-#include <dumux/discretization/box/fvelementgeometryvector.hh>
+#include <dumux/discretization/box/globalfvgeometry.hh>
+#include <dumux/discretization/box/fvelementgeometry.hh>
#include <dumux/discretization/box/stencils.hh>
#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
@@ -58,7 +58,10 @@ namespace Properties {
SET_INT_PROP(BoxModel, DiscretizationMethod, GET_PROP(TypeTag, DiscretizationMethods)::Box);
//! Set the default for the FVElementGeometry vector
-SET_TYPE_PROP(BoxModel, FVElementGeometryVector, BoxFVElementGeometryVector<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVElementGeometryCache)>);
+SET_TYPE_PROP(BoxModel, GlobalFVGeometry, BoxGlobalFVGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
+
+//! Set the default for the FVElementGeometry vector
+SET_TYPE_PROP(BoxModel, FVElementGeometry, BoxFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
//! The sub control volume
SET_PROP(BoxModel, SubControlVolume)
@@ -71,7 +74,7 @@ private:
using ScvGeometry = Dune::MultiLinearGeometry<Scalar, dim, dimWorld>;
using IndexType = typename GridView::IndexSet::IndexType;
public:
- using type = BoxSubControlVolume<ScvGeometry, IndexType, /*isBox=*/true>;
+ using type = BoxSubControlVolume<ScvGeometry, IndexType>;
};
SET_PROP(BoxModel, SubControlVolumeFace)
@@ -103,7 +106,7 @@ SET_TYPE_PROP(BoxModel, CurrentVolumeVariablesVector, BoxVolumeVariablesVector<T
SET_TYPE_PROP(BoxModel, PreviousVolumeVariablesVector, BoxVolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
//! The global flux variables cache vector class
-SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, BoxFluxVariablesCacheVector<TypeTag>);
+SET_TYPE_PROP(BoxModel, FluxVariablesCacheVector, BoxFluxVariablesCacheVector<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
//! Set the BaseLocalResidual to BoxLocalResidual
SET_TYPE_PROP(BoxModel, BaseLocalResidual, BoxLocalResidual<TypeTag>);
diff --git a/dumux/implicit/cellcentered/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh
index a11630ed96..addff4924c 100644
--- a/dumux/implicit/cellcentered/elementboundarytypes.hh
+++ b/dumux/implicit/cellcentered/elementboundarytypes.hh
@@ -87,12 +87,12 @@ public:
(*this)[0].reset();
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
if (!problem.model().onBoundary(scv))
return;
- for (const auto& scvFace : scvfs(fvGeometry))
+ for (auto&& scvFace : scvfs(fvGeometry))
{
if (!scvFace.boundary())
continue;
diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh
index a38fe2d640..3ddce0a1cb 100644
--- a/dumux/implicit/cellcentered/localjacobian.hh
+++ b/dumux/implicit/cellcentered/localjacobian.hh
@@ -75,7 +75,7 @@ class CCLocalJacobian : public ImplicitLocalJacobian<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
@@ -114,8 +114,8 @@ public:
assemblyMap_.resize(problem.gridView().size(0));
for (const auto& element : elements(problem.gridView()))
{
- // bind FVGeometry to the element to prepare all the geometries of the stencil
- problem.model().fvGeometries_().bind(element);
+ // get a local finite volume geometry object that is bindable
+ auto fvGeometryJ = localView(problem.model().globalFvGeometry());
auto globalI = problem.elementMapper().index(element);
const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
@@ -123,18 +123,19 @@ public:
assemblyMap_[globalI].reserve(neighborStencil.size());
for (auto globalJ : neighborStencil)
{
- const auto& elementJ = this->model_().fvGeometries().element(globalJ);
- const auto& fvGeometry = this->model_().fvGeometries(elementJ);
+ const auto& elementJ = fvGeometryJ.globalFvGeometry().element(globalJ);
+ fvGeometryJ.bind(elementJ);
// find the flux vars needed for the calculation of the flux into element
std::vector<IndexType> fluxVarIndices;
- for (const auto& scvFaceJ : scvfs(fvGeometry))
+ for (auto&& scvFaceJ : scvfs(fvGeometryJ))
{
auto fluxVarsIdx = scvFaceJ.index();
// if globalI is in flux var stencil, add to list
FluxVariables fluxVars;
- const auto fluxStencil = fluxVars.computeStencil(problem, elementJ, scvFaceJ);
+ const auto fluxStencil = fluxVars.computeStencil(problem, elementJ, fvGeometryJ, scvFaceJ);
+
for (auto globalIdx : fluxStencil)
if (globalIdx == globalI)
fluxVarIndices.push_back(fluxVarsIdx);
@@ -156,15 +157,16 @@ public:
const bool isGhost = (element.partitionType() == Dune::GhostEntity);
// prepare the volvars/fvGeometries in case caching is disabled
- this->model_().fvGeometries_().bind(element);
- this->model_().curVolVars_().bind(element);
- this->model_().prevVolVars_().bindElement(element);
- this->model_().fluxVariablesCache_().bind(element);
+ auto fvGeometry = localView(this->model_().globalFvGeometry());
+ fvGeometry.bind(element);
+
+ this->model_().curVolVars_().bind(element, fvGeometry);
+ this->model_().prevVolVars_().bindElement(element, fvGeometry);
+ this->model_().fluxVariablesCache_().bind(element, fvGeometry);
// set the actual dof index
globalI_ = this->problem_().elementMapper().index(element);
- const auto& fvGeometry = this->model_().fvGeometries(element);
bcTypes_.update(this->problem_(), element, fvGeometry);
// calculate the local residual
@@ -175,7 +177,7 @@ public:
}
else
{
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes_, fvGeometry);
this->residual_ = this->localResidual().residual();
// store residual in global container as well
residual[globalI_] = this->localResidual().residual(0);
@@ -184,7 +186,7 @@ public:
this->model_().updatePVWeights(fvGeometry);
// calculate derivatives of all dofs in stencil with respect to the dofs in the element
- evalPartialDerivatives_(element, matrix, residual, isGhost);
+ evalPartialDerivatives_(element, fvGeometry, matrix, residual, isGhost);
// TODO: calculate derivatives in the case of an extended source stencil
// const auto& extendedSourceStencil = model_().stencils(element).extendedSourceStencil();
@@ -205,7 +207,11 @@ public:
{ return assemblyMap_; }
private:
- void evalPartialDerivatives_(const Element& element, JacobianMatrix& matrix, SolutionVector& residual, const bool isGhost)
+ void evalPartialDerivatives_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ JacobianMatrix& matrix,
+ SolutionVector& residual,
+ const bool isGhost)
{
// get stencil informations
const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
@@ -230,16 +236,19 @@ private:
unsigned int j = 0;
for (auto globalJ : neighborStencil)
{
- auto elementJ = this->problem_().model().fvGeometries().element(globalJ);
+ auto elementJ = fvGeometry.globalFvGeometry().element(globalJ);
neighborElements.push_back(elementJ);
for (auto fluxVarIdx : assemblyMap_[globalI_][j])
- origFlux[j] += localRes.evalFlux_(elementJ, fluxVarIdx);
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, scvf);
+ }
++j;
}
- const auto& scv = this->model_().fvGeometries().subControlVolume(globalI_);
+ auto&& scv = fvGeometry.scv(globalI_);
VolumeVariables origVolVars(this->model_().curVolVars(scv));
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
@@ -261,12 +270,12 @@ private:
// update the volume variables and bind the flux var cache again
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
- this->model_().fluxVariablesCache_().bind(element);
+ this->model_().fluxVariablesCache_().bind(element, fvGeometry);
if (!isGhost)
{
// calculate the residual with the deflected primary variables
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes_, fvGeometry);
// store the residual and the storage term
partialDeriv = this->localResidual().residual(0);
@@ -276,7 +285,10 @@ private:
for (std::size_t k = 0; k < numNeighbors; ++k)
{
for (auto fluxVarIdx : assemblyMap_[globalI_][k])
- neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fluxVarIdx);
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, scvf);
+ }
}
}
else
@@ -300,12 +312,12 @@ private:
// update the volume variables and bind the flux var cache again
this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
- this->model_().fluxVariablesCache_().bind(element);
+ this->model_().fluxVariablesCache_().bind(element, fvGeometry);
if (!isGhost)
{
// calculate the residual with the deflected primary variables
- this->localResidual().eval(element, bcTypes_);
+ this->localResidual().eval(element, bcTypes_, fvGeometry);
// subtract the residual from the derivative storage
partialDeriv -= this->localResidual().residual(0);
@@ -315,7 +327,10 @@ private:
for (std::size_t k = 0; k < numNeighbors; ++k)
{
for (auto fluxVarIdx : assemblyMap_[globalI_][k])
- neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fluxVarIdx);
+ {
+ auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+ neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, scvf);
+ }
}
}
else
diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh
index 4766debb5f..2882022e1d 100644
--- a/dumux/implicit/cellcentered/localresidual.hh
+++ b/dumux/implicit/cellcentered/localresidual.hh
@@ -43,24 +43,21 @@ namespace Dumux
template<class TypeTag>
class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
{
- typedef ImplicitLocalResidual<TypeTag> ParentType;
+ using ParentType = ImplicitLocalResidual<TypeTag>;
friend class ImplicitLocalResidual<TypeTag>;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, NumEq)
- };
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) };
-
enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache) };
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
// copying the local residual class is not a good idea
@@ -70,42 +67,46 @@ public:
protected:
- void evalFluxes_()
+ void evalFluxes_(const Element& element,
+ const ElementBoundaryTypes& bcTypes,
+ const FVElementGeometry& fvGeometry)
{
// calculate the mass flux over the scv faces
- const auto& fvGeometry = this->fvGeometry_();
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
- this->residual_[0] += this->asImp_().computeFlux_(scvf);
+ this->residual_[0] += this->asImp_().computeFlux_(element, fvGeometry, scvf, bcTypes);
}
}
- PrimaryVariables computeFlux_(const SubControlVolumeFace &scvf)
+ PrimaryVariables computeFlux_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace &scvf,
+ const ElementBoundaryTypes& bcTypes)
{
if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
- return this->asImp_().computeFlux(scvf);
+ return this->asImp_().computeFlux(element, fvGeometry, scvf, bcTypes[0]);
else
return PrimaryVariables(0.0);
}
- void evalBoundary_()
+ void evalBoundary_(const Element &element,
+ const ElementBoundaryTypes& bcTypes,
+ const FVElementGeometry& fvGeometry)
{
- const auto& fvGeometry = this->fvGeometry_();
-
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
- this->residual_[0] += evalBoundaryFluxes_(scvf);
+ this->residual_[0] += evalBoundaryFluxes_(element, scvf, fvGeometry, bcTypes[0]);
}
// additionally treat mixed D/N conditions in a strong sense
- if (this->bcTypes_().hasDirichlet())
+ if (bcTypes.hasDirichlet())
{
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
- this->asImp_().evalDirichlet_(scvf);
+ this->asImp_().evalDirichlet_(element, scvf, fvGeometry);
}
}
}
@@ -114,14 +115,15 @@ protected:
* \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet
* boundary conditions to the local residual.
*/
- PrimaryVariables evalBoundaryFluxes_(const SubControlVolumeFace &scvf)
+ PrimaryVariables evalBoundaryFluxes_(const Element &element,
+ const SubControlVolumeFace &scvf,
+ const FVElementGeometry& fvGeometry,
+ const BoundaryTypes& bcTypes)
{
- BoundaryTypes bcTypes = this->problem_().boundaryTypes(this->element_(), scvf);
-
// evaluate the Neumann conditions at the boundary face
PrimaryVariables flux(0);
if (bcTypes.hasNeumann() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
- flux += this->asImp_().evalNeumannSegment_(scvf, bcTypes);
+ flux += this->asImp_().evalNeumannSegment_(element, scvf, fvGeometry, bcTypes);
// TODO: evaluate the outflow conditions at the boundary face
//if (bcTypes.hasOutflow() /*TODO: && !GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
@@ -129,37 +131,42 @@ protected:
// evaluate the pure Dirichlet conditions at the boundary face
if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
- flux += this->asImp_().evalDirichletSegment_(scvf, bcTypes);
+ flux += this->asImp_().evalDirichletSegment_(element, scvf, fvGeometry, bcTypes);
return flux;
}
+
/*!
* \brief Evaluate Dirichlet conditions on faces that have mixed
* Dirichlet/Neumann boundary conditions.
*/
- void evalDirichlet_(const SubControlVolumeFace &scvf)
+ void evalDirichlet_(const Element& element,
+ const SubControlVolumeFace &scvf,
+ const FVElementGeometry& fvGeometry)
{
- BoundaryTypes bcTypes = this->problem_().boundaryTypes(this->element_(), scvf);
+ BoundaryTypes bcTypes = this->problem().boundaryTypes(element, scvf);
if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
- this->asImp_().evalDirichletSegmentMixed_(scvf, bcTypes);
+ this->asImp_().evalDirichletSegmentMixed_(element, scvf, fvGeometry, bcTypes);
}
/*!
* \brief Add Neumann boundary conditions for a single scv face
*/
- PrimaryVariables evalNeumannSegment_(const SubControlVolumeFace &scvf,
+ PrimaryVariables evalNeumannSegment_(const Element& element,
+ const SubControlVolumeFace &scvf,
+ const FVElementGeometry& fvGeometry,
const BoundaryTypes &bcTypes)
{
// temporary vector to store the neumann boundary fluxes
PrimaryVariables flux(0);
- auto neumannFluxes = this->problem_().neumann(this->element_(), scvf);
+ auto neumannFluxes = this->problem().neumann(element, scvf);
// multiply neumann fluxes with the area and the extrusion factor
- const auto& scv = this->problem_().model().fvGeometries().subControlVolume(scvf.insideScvIdx());
- neumannFluxes *= scvf.area()*this->problem_().model().curVolVars(scv).extrusionFactor();
+ auto&& scv = fvGeometry.scv(scvf.insideScvIdx());
+ neumannFluxes *= scvf.area()*this->problem().model().curVolVars(scv).extrusionFactor();
// add fluxes to the residual
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@@ -169,6 +176,66 @@ protected:
return flux;
}
+ /*!
+ * \brief Treat Dirichlet boundary conditions in a weak sense for a single
+ * intersection that only has Dirichlet boundary conditions
+ */
+ PrimaryVariables evalDirichletSegment_(const Element& element,
+ const SubControlVolumeFace &scvf,
+ const FVElementGeometry& fvGeometry,
+ const BoundaryTypes &bcTypes)
+ {
+ // temporary vector to store the Dirichlet boundary fluxes
+ PrimaryVariables flux(0);
+
+ if (!constantBC || !enableVolVarsCache)
+ {
+ // update corresponding boundary volume variables before flux calculation
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+ const auto dirichletPriVars = this->problem().dirichlet(element, scvf);
+ this->problem().model().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars,
+ this->problem(),
+ element,
+ insideScv);
+ }
+
+ auto dirichletFlux = this->asImp_().computeFlux(element, fvGeometry, scvf, bcTypes);
+
+ // add fluxes to the residual
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ if (bcTypes.isDirichlet(eqIdx))
+ flux[eqIdx] += dirichletFlux[eqIdx];
+
+ return flux;
+ }
+
+ /*!
+ * \brief Treat Dirichlet boundary conditions in a strong sense for a
+ * single intersection that has mixed D/N boundary conditions
+ */
+ void evalDirichletSegmentMixed_(const Element& element,
+ const SubControlVolumeFace &scvf,
+ const FVElementGeometry& fvGeometry,
+ const BoundaryTypes &bcTypes)
+ {
+ // temporary vector to store the Dirichlet values
+ PrimaryVariables dirichletValues = this->problem().dirichlet(element, scvf);
+
+ // get the primary variables
+ const auto& priVars = this->problem().model().curVolVars(scvf.insideScvIdx()).priVars();
+
+ // set Dirichlet conditions in a strong sense
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isDirichlet(eqIdx))
+ {
+ auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+ }
+ }
+ }
+
/*!
* \brief Add outflow boundary conditions for a single intersection
*/
@@ -227,59 +294,6 @@ protected:
// restore the pointer to the FVElementGeometry
this->fvElemGeomPtr_ = oldFVGeometryPtr;
}*/
-
- /*!
- * \brief Treat Dirichlet boundary conditions in a weak sense for a single
- * intersection that only has Dirichlet boundary conditions
- */
- PrimaryVariables evalDirichletSegment_(const SubControlVolumeFace &scvf,
- const BoundaryTypes &bcTypes)
- {
- // temporary vector to store the Dirichlet boundary fluxes
- PrimaryVariables flux(0);
-
- if (!constantBC || !enableVolVarsCache)
- {
- // update corresponding boundary volume variables before flux calculation
- const auto insideScvIdx = scvf.insideScvIdx();
- const auto& insideScv = this->model_().fvGeometries().subControlVolume(insideScvIdx);
- const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf);
- this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv);
- }
-
- auto dirichletFlux = this->asImp_().computeFlux(scvf);
-
- // add fluxes to the residual
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- if (bcTypes.isDirichlet(eqIdx))
- flux[eqIdx] += dirichletFlux[eqIdx];
-
- return flux;
- }
-
- /*!
- * \brief Treat Dirichlet boundary conditions in a strong sense for a
- * single intersection that has mixed D/N boundary conditions
- */
- void evalDirichletSegmentMixed_(const SubControlVolumeFace &scvf,
- const BoundaryTypes &bcTypes)
- {
- // temporary vector to store the Dirichlet values
- PrimaryVariables dirichletValues = this->problem_().dirichlet(this->element_(), scvf);
-
- // get the primary variables
- const auto& priVars = this->model_().curVolVars(scvf.insideScvIdx()).priVars();
-
- // set Dirichlet conditions in a strong sense
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isDirichlet(eqIdx))
- {
- auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- this->residual_[0][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
- }
};
}
diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
index cc904e2b56..5ce873a384 100644
--- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
@@ -29,7 +29,8 @@
#include <dumux/implicit/propertydefaults.hh>
#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
-#include <dumux/discretization/cellcentered/tpfa/fvelementgeometryvector.hh>
+#include <dumux/discretization/cellcentered/tpfa/globalfvgeometry.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
#include <dumux/discretization/cellcentered/tpfa/subcontrolvolume.hh>
#include <dumux/discretization/cellcentered/tpfa/subcontrolvolumeface.hh>
#include <dumux/implicit/cellcentered/properties.hh>
@@ -43,8 +44,11 @@ namespace Properties {
//! Set the corresponding discretization method property
SET_INT_PROP(CCTpfaModel, DiscretizationMethod, GET_PROP(TypeTag, DiscretizationMethods)::CCTpfa);
-//! Set the default for the FVElementGeometry vector
-SET_TYPE_PROP(CCTpfaModel, FVElementGeometryVector, CCTpfaFVElementGeometryVector<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVElementGeometryCache)>);
+//! Set the default for the global finite volume geometry
+SET_TYPE_PROP(CCTpfaModel, GlobalFVGeometry, CCTpfaGlobalFVGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
+
+//! Set the default for the local finite volume geometry
+SET_TYPE_PROP(CCTpfaModel, FVElementGeometry, CCTpfaFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
//! The sub control volume
SET_PROP(CCTpfaModel, SubControlVolume)
diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh
index 96c8ee5cc3..2396859acd 100644
--- a/dumux/implicit/localresidual.hh
+++ b/dumux/implicit/localresidual.hh
@@ -51,7 +51,7 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
@@ -82,6 +82,14 @@ public:
void init(Problem &problem)
{ problemPtr_ = &problem; }
+
+ /*!
+ * \name User interface
+ * \note The following methods are usually expensive to evaluate
+ * They are useful for outputting residual information.
+ */
+ // \{
+
/*!
* \brief Compute the local residual, i.e. the deviation of the
* equations from zero.
@@ -92,15 +100,17 @@ public:
void eval(const Element &element)
{
// make sure FVElementGeometry and volume variables are bound to the element
- model_().fvGeometries_().bind(element);
- model_().curVolVars_().bind(element);
- model_().prevVolVars_().bindElement(element);
- model_().fluxVariablesCache_().bindElement(element);
+ auto fvGeometry = localView(problem().model().globalFvGeometry());
+ fvGeometry.bind(element);
+
+ problem().model().curVolVars_().bind(element, fvGeometry);
+ problem().model().prevVolVars_().bindElement(element, fvGeometry);
+ problem().model().fluxVariablesCache_().bindElement(element, fvGeometry);
ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
+ bcTypes.update(problem(), element, fvGeometry);
- asImp_().eval(element, bcTypes);
+ asImp_().eval(element, bcTypes, fvGeometry);
}
/*!
@@ -114,47 +124,53 @@ public:
*/
void evalStorage(const Element &element)
{
- elemPtr_ = &element;
-
// make sure FVElementGeometry and volume variables are bound to the element
- model_().fvGeometries_().bindElement(element);
- model_().curVolVars_().bindElement(element);
- model_().prevVolVars_().bindElement(element);
+ auto fvGeometry = localView(problem().model().globalFvGeometry());
+ fvGeometry.bind(element);
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
+ // make sure FVElementGeometry and volume variables are bound to the element
+ problem().model().curVolVars_().bindElement(element, fvGeometry);
+ problem().model().prevVolVars_().bindElement(element, fvGeometry);
- asImp_().evalStorage_();
+ asImp_().evalStorage_(fvGeometry);
}
- /*!
- * \brief Compute the flux term for the current solution.
- *
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param curVolVars The volume averaged variables for all
- * sub-contol volumes of the element
- */
- void evalFluxes(const Element &element)
- {
- elemPtr_ = &element;
+ // !
+ // * \brief Compute the flux term for the current solution.
+ // *
+ // * \param element The DUNE Codim<0> entity for which the residual
+ // * ought to be calculated
+ // * \param curVolVars The volume averaged variables for all
+ // * sub-contol volumes of the element
- // make sure FVElementGeometry and volume variables are bound to the element
- model_().fvGeometries_().bind(element);
- model_().curVolVars_().bind(element);
- model_().prevVolVars_().bindElement(element);
- model_().fluxVariablesCache_().bindElement(element);
+ // void evalFluxes(const Element &element)
+ // {
+ // elemPtr_ = &element;
- ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
+ // // make sure FVElementGeometry and volume variables are bound to the element
+ // problem().model().fvGeometries_().bind(element);
+ // problem().model().curVolVars_().bind(element);
+ // problem().model().prevVolVars_().bindElement(element);
+ // problem().model().fluxVariablesCache_().bindElement(element);
- residual_.resize(fvGeometry_().numScv);
- residual_ = 0;
+ // ElementBoundaryTypes bcTypes;
+ // bcTypes.update(problem(), element, fvGeometry_());
+
+ // residual_.resize(fvGeometry_().numScv);
+ // residual_ = 0;
+
+ // bcTypesPtr_ = &bcTypes;
+ // asImp_().evalFluxes_();
+ // }
+
+ // \}
- bcTypesPtr_ = &bcTypes;
- asImp_().evalFluxes_();
- }
+
+ /*!
+ * \name Main interface
+ * \note Methods used by the assembler to compute derivatives and residual
+ */
+ // \{
/*!
* \brief Compute the local residual, i.e. the deviation of the
@@ -173,40 +189,38 @@ public:
* vertices of the element
*/
void eval(const Element &element,
- const ElementBoundaryTypes &bcTypes)
+ const ElementBoundaryTypes &bcTypes,
+ const FVElementGeometry& fvGeometry)
{
- // At this point the fv geometry has to be bound already
- elemPtr_ = &element;
- bcTypesPtr_ = &bcTypes;
-
// resize the vectors for all terms
- int numScv = fvGeometry_().numScv();
+ auto numScv = fvGeometry.numScv();
residual_.resize(numScv);
storageTerm_.resize(numScv);
residual_ = 0.0;
storageTerm_ = 0.0;
- asImp_().evalFluxes_();
- asImp_().evalVolumeTerms_();
- asImp_().evalBoundary_();
+ asImp_().evalFluxes_(element, bcTypes, fvGeometry);
+ asImp_().evalVolumeTerms_(element, bcTypes, fvGeometry);
+ asImp_().evalBoundary_(element, bcTypes, fvGeometry);
}
/*!
* \brief Calculate the source term of the equation
*
- * \param source The source/sink in the sub-control volume for each phase
- * \param scvIdx The index of the sub-control volume
+ * \param scv The sub-control volume over which we integrate the source term
*
*/
- PrimaryVariables computeSource(const SubControlVolume &scv)
+ PrimaryVariables computeSource(const Element& element,
+ const SubControlVolume &scv)
{
PrimaryVariables source(0);
- source += this->problem_().source(element_(), scv);
+ // add contributions from volume flux sources
+ source += this->problem().source(element, scv);
// add contribution from possible point sources
- source += this->problem_().scvPointSources(element_(), scv);
+ source += this->problem().scvPointSources(element, scv);
return source;
}
@@ -241,80 +255,75 @@ public:
const PrimaryVariables &storageTerm(const int scvIdx) const
{ return storageTerm_[scvIdx]; }
+ /*!
+ * \brief Return the problem we are solving. Only call this after init()!
+ */
+ const Problem& problem() const
+ { return *problemPtr_; }
+
+ /*!
+ * \brief Return the problem we are solving. Only call this after init()!
+ */
+ Problem& problem()
+ { return *problemPtr_; }
+
+ // \}
+
protected:
Implementation &asImp_()
- {
- assert(static_cast<Implementation*>(this) != 0);
- return *static_cast<Implementation*>(this);
- }
+ { return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
- {
- assert(static_cast<const Implementation*>(this) != 0);
- return *static_cast<const Implementation*>(this);
- }
+ { return *static_cast<const Implementation*>(this); }
- PrimaryVariables evalFlux_(const Element &element, const int scvfIdx)
+ PrimaryVariables evalFlux_(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace scvf)
{
- const auto& scvf = model_().fvGeometries().subControlVolumeFace(scvfIdx);
- elemPtr_ = &element;
-
ElementBoundaryTypes bcTypes;
- bcTypes.update(problem_(), element, fvGeometry_());
- bcTypesPtr_ = &bcTypes;
+ bcTypes.update(problem(), element, fvGeometry);
- residual_.resize(fvGeometry_().numScv());
+ residual_.resize(fvGeometry.numScv());
residual_ = 0;
- return evalFlux_(scvf);
- }
-
- PrimaryVariables evalFlux_(const int scvfIdx)
- {
- auto&& scvf = model_().fvGeometries().subControlVolumeFace(scvfIdx);
- return evalFlux_(scvf);
- }
-
- PrimaryVariables evalFlux_(const SubControlVolumeFace &scvf)
- {
- return asImp_().computeFlux_(scvf);
+ return asImp_().computeFlux_(element, fvGeometry, scvf, bcTypes);
}
/*!
* \brief Set the local residual to the storage terms of all
* sub-control volumes of the current element.
*/
- void evalStorage_()
+ void evalStorage_(const FVElementGeometry& fvGeometry)
{
- storageTerm_.resize(fvGeometry_().numScv());
+ storageTerm_.resize(fvGeometry.numScv());
storageTerm_ = 0;
// calculate the amount of conservation each quantity inside
// all sub control volumes
- const auto& fvGeometry = fvGeometry_();
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
- auto scvIdx = scv.indexInElement();
+ auto scvIdx = isBox ? scv.index() : 0;
- storageTerm_[scvIdx] = asImp_().computeStorage(scv, model_().curVolVars(scv));
- storageTerm_[scvIdx] *= scv.volume() * model_().curVolVars(scv).extrusionFactor();
+ storageTerm_[scvIdx] = asImp_().computeStorage(scv, problem().model().curVolVars(scv));
+ storageTerm_[scvIdx] *= scv.volume() * problem().model().curVolVars(scv).extrusionFactor();
}
}
- PrimaryVariables evalSource_()
- {
- PrimaryVariables source(0);
- const auto& fvGeometry = fvGeometry_();
- for (const auto& scv : scvs(fvGeometry))
- {
- source += this->problem_().source(element_(), scv);
+ // PrimaryVariables evalSource_(const Element& element,
+ // const FVElementGeometry& fvGeometry)
+ // {
+ // PrimaryVariables source(0);
+ // const auto& fvGeometry = fvGeometry.fvElementGeometry();
+ // for (auto&& scv : scvs(fvGeometry))
+ // {
+ // source += this->problem().source(element, scv);
- // add contribution from possible point sources
- source += this->problem_().scvPointSources(element_(), scv);
- }
+ // // add contribution from possible point sources
+ // source += this->problem().scvPointSources(element, scv);
+ // }
- return source;
- }
+ // return source;
+ // }
/*!
* \brief Add the change the source term for stationary problems
@@ -323,17 +332,18 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTerms_()
+ evalVolumeTerms_(const Element &element,
+ const ElementBoundaryTypes &bcTypes,
+ const FVElementGeometry& fvGeometry)
{
// evaluate the volume terms (storage + source terms)
- const auto& fvGeometry = fvGeometry_();
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
- auto scvIdx = scv.indexInElement();
- auto curExtrusionFactor = model_().curVolVars(scv).extrusionFactor();
+ auto scvIdx = isBox ? scv.index() : 0;
+ auto curExtrusionFactor = problem().model().curVolVars(scv).extrusionFactor();
// subtract the source term from the local rate
- PrimaryVariables source = asImp_().computeSource(scv);
+ PrimaryVariables source = asImp_().computeSource(element, scv);
source *= scv.volume()*curExtrusionFactor;
residual_[scvIdx] -= source;
@@ -347,15 +357,16 @@ protected:
*/
template<class P = Problem>
typename std::enable_if<!Dumux::Capabilities::isStationary<P>::value, void>::type
- evalVolumeTerms_()
+ evalVolumeTerms_(const Element &element,
+ const ElementBoundaryTypes &bcTypes,
+ const FVElementGeometry& fvGeometry)
{
// evaluate the volume terms (storage + source terms)
- const auto& fvGeometry = fvGeometry_();
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
- auto scvIdx = scv.indexInElement();
- auto prevExtrusionFactor = model_().prevVolVars(scv).extrusionFactor();
- auto curExtrusionFactor = model_().curVolVars(scv).extrusionFactor();
+ auto scvIdx = isBox ? scv.index() : 0;
+ auto prevExtrusionFactor = problem().model().prevVolVars(scv).extrusionFactor();
+ auto curExtrusionFactor = problem().model().curVolVars(scv).extrusionFactor();
// mass balance within the element. this is the
// \f$\frac{m}{\partial t}\f$ term if using implicit
@@ -363,8 +374,8 @@ protected:
//
// We might need a more explicit way for
// doing the time discretization...
- PrimaryVariables prevStorage = asImp_().computeStorage(scv, model_().prevVolVars(scv));
- PrimaryVariables curStorage = asImp_().computeStorage(scv, model_().curVolVars(scv));
+ PrimaryVariables prevStorage = asImp_().computeStorage(scv, problem().model().prevVolVars(scv));
+ PrimaryVariables curStorage = asImp_().computeStorage(scv, problem().model().curVolVars(scv));
prevStorage *= prevExtrusionFactor;
curStorage *= curExtrusionFactor;
@@ -372,96 +383,25 @@ protected:
storageTerm_[scvIdx] = std::move(curStorage);
storageTerm_[scvIdx] -= std::move(prevStorage);
storageTerm_[scvIdx] *= scv.volume();
- storageTerm_[scvIdx] /= problem_().timeManager().timeStepSize();
+ storageTerm_[scvIdx] /= problem().timeManager().timeStepSize();
// add the storage term to the residual
residual_[scvIdx] += storageTerm_[scvIdx];
// subtract the source term from the local rate
- PrimaryVariables source = asImp_().computeSource(scv);
+ PrimaryVariables source = asImp_().computeSource(element, scv);
source *= scv.volume()*curExtrusionFactor;
residual_[scvIdx] -= source;
}
}
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem_() const
- { return *problemPtr_; }
-
- /*!
- * \brief Returns a reference to the problem.
- */
- Problem &problem_()
- { return *problemPtr_; }
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model_() const
- { return problem_().model(); }
-
- /*!
- * \brief Returns a reference to the model.
- */
- Model &model_()
- { return problem_().model(); }
-
- /*!
- * \brief Returns a reference to the grid view.
- */
- const GridView &gridView_() const
- { return problem_().gridView(); }
-
- /*!
- * \brief Returns a reference to the current element.
- */
- const Element &element_() const
- {
- Valgrind::CheckDefined(elemPtr_);
- return *elemPtr_;
- }
-
- /*!
- * \brief Returns a reference to the current element's finite
- * volume geometry.
- */
- const FVElementGeometry& fvGeometry_() const
- {
- return model_().fvGeometries(element_());
- }
-
- /*!
- * \brief Returns a reference to the boundary types of all
- * sub-control volumes of the current element.
- */
- const ElementBoundaryTypes &bcTypes_() const
- {
- Valgrind::CheckDefined(bcTypesPtr_);
- return *bcTypesPtr_;
- }
-
- /*!
- * \brief Returns a reference to the boundary types of the i-th
- * sub-control volume of the current element.
- */
- const BoundaryTypes &bcTypes_(const int i) const
- {
- return bcTypes_()[i];
- }
-
protected:
ElementSolutionVector storageTerm_;
ElementSolutionVector residual_;
- // The problem we would like to solve
- Problem *problemPtr_;
-
- const Element *elemPtr_;
-
- const ElementBoundaryTypes *bcTypesPtr_;
+private:
+ Problem* problemPtr_;
};
}
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index 8a76072f71..4088bf1cc9 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -82,7 +82,7 @@ class ImplicitModel
dim = GridView::dimension
};
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometryVector) FVElementGeometryVector;
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) LocalJacobian;
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
@@ -119,8 +119,8 @@ public:
updateBoundaryIndices_();
- fvGeometryVector_ = std::make_shared<FVElementGeometryVector>(gridView_());
- fvGeometryVector_->update(problem_());
+ globalFvGeometryPtr_ = std::make_shared<GlobalFVGeometry>(gridView_());
+ globalFvGeometryPtr_->update(problem);
int numDofs = asImp_().numDofs();
uCur_.resize(numDofs);
@@ -132,18 +132,18 @@ public:
asImp_().applyInitialSolution_();
// resize and update the volVars with the initial solution
- curVolVarsVector_.update(problem_(), curSol());
+ curVolVarsVector_.update(problem, curSol());
// update stencils
- stencilsVector_.update(problem_());
+ stencilsVector_.update(problem);
// update the flux variables caches
- fluxVarsCacheVector_.update(problem_());
+ fluxVarsCacheVector_.update(problem);
// initialize assembler and create matrix
- localJacobian_.init(problem_());
+ localJacobian_.init(problem);
jacAsm_ = std::make_shared<JacobianAssembler>();
- jacAsm_->init(problem_());
+ jacAsm_->init(problem);
// also set the solution of the "previous" time step to the
// initial solution.
@@ -786,11 +786,8 @@ public:
const VolumeVariables& curVolVars(const unsigned int scvIdx) const
{ return curVolVarsVector_[scvIdx]; }
- const FVElementGeometryVector& fvGeometries() const
- { return *fvGeometryVector_; }
-
- const FVElementGeometry& fvGeometries(const Element& element) const
- { return fvGeometryVector_->fvGeometry(element); }
+ const GlobalFVGeometry& globalFvGeometry() const
+ { return *globalFvGeometryPtr_; }
protected:
@@ -815,8 +812,8 @@ protected:
FluxVariablesCache& fluxVarsCache_(const SubControlVolumeFace& scvf)
{ return fluxVarsCacheVector_[scvf]; }
- FVElementGeometryVector& fvGeometries_()
- { return *fvGeometryVector_; }
+ GlobalFVGeometry& globalFvGeometry_()
+ { return *globalFvGeometryPtr_; }
/*!
* \brief A reference to the problem on which the model is applied.
@@ -858,11 +855,11 @@ protected:
for (const auto& element : elements(gridView_()))
{
// deal with the current element, bind FVGeometry to it
- fvGeometries_().bindElement(element);
- const auto& fvGeometry = fvGeometries(element);
+ auto fvGeometry = localView(globalFvGeometry());
+ fvGeometry.bindElement(element);
// loop over sub control volumes
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
// let the problem do the dirty work of nailing down
// the initial solution.
@@ -980,7 +977,7 @@ protected:
StencilsVector stencilsVector_;
// the finite volume element geometries
- std::shared_ptr<FVElementGeometryVector> fvGeometryVector_;
+ std::shared_ptr<GlobalFVGeometry> globalFvGeometryPtr_;
// container to store the box volumes
Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
diff --git a/dumux/implicit/problem.hh b/dumux/implicit/problem.hh
index 77c9cd2987..916cb385dd 100644
--- a/dumux/implicit/problem.hh
+++ b/dumux/implicit/problem.hh
@@ -973,8 +973,8 @@ public:
PrimaryVariables scvPointSources(const Element &element, const SubControlVolume &scv) const
{
PrimaryVariables source(0);
-
- auto key = std::make_pair(this->gridView().indexSet().index(element), scv.indexInElement());
+ auto scvIdx = isBox ? scv.index() : 0;
+ auto key = std::make_pair(this->gridView().indexSet().index(element), scvIdx);
if (pointSourceMap_.count(key))
{
// call the solDependent function. Herein the user might fill/add values to the point sources
diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh
index c4effc39f9..cff3fdf63c 100644
--- a/dumux/implicit/properties.hh
+++ b/dumux/implicit/properties.hh
@@ -59,9 +59,9 @@ NEW_PROP_TAG(LocalJacobian); //!< The type of the local jacobian operator
NEW_PROP_TAG(SubControlVolume);//!< The type of the sub control volume
NEW_PROP_TAG(SubControlVolumeFace); //!< The type of the sub control volume face
-NEW_PROP_TAG(FVElementGeometry); //!< The type of the finite volume geometry (iterators over scvs, scvfs)
-NEW_PROP_TAG(FVElementGeometryVector); //!< The type of the finite volume geometry vector
-NEW_PROP_TAG(EnableGlobalFVElementGeometryCache); //! specifies if geometric data is be saved (faster, but more memory consuming)
+NEW_PROP_TAG(FVElementGeometry); //!< The type of the local finite volume geometry (iterators over scvs, scvfs)
+NEW_PROP_TAG(GlobalFVGeometry); //!< The type of the global finite volume geometry
+NEW_PROP_TAG(EnableGlobalFVGeometryCache); //! specifies if geometric data is be saved (faster, but more memory consuming)
NEW_PROP_TAG(JacobianAssembler); //!< Assembles the global jacobian matrix
NEW_PROP_TAG(JacobianMatrix); //!< Type of the global jacobian matrix
diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh
index 2d4219e401..864a2b0202 100644
--- a/dumux/implicit/propertydefaults.hh
+++ b/dumux/implicit/propertydefaults.hh
@@ -39,7 +39,6 @@
#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
#include <dumux/discretization/volumevariables.hh>
-#include <dumux/discretization/fvelementgeometry.hh>
#include <dumux/discretization/darcyslaw.hh>
#include <dumux/discretization/fickslaw.hh>
@@ -106,9 +105,6 @@ SET_TYPE_PROP(ImplicitBase,
//! Set the BaseModel to ImplicitModel
SET_TYPE_PROP(ImplicitBase, BaseModel, ImplicitModel<TypeTag>);
-//! The finite volume element geometry providing iterators over scvs and scv faces
-SET_TYPE_PROP(ImplicitBase, FVElementGeometry, Dumux::FVElementGeometry<TypeTag>);
-
//! The volume variable class, to be overloaded by the model
SET_TYPE_PROP(ImplicitBase, VolumeVariables, ImplicitVolumeVariables<TypeTag>);
@@ -167,7 +163,7 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false);
SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false);
//! We do not store the FVGeometry by default
-SET_BOOL_PROP(ImplicitBase, EnableGlobalFVElementGeometryCache, false);
+SET_BOOL_PROP(ImplicitBase, EnableGlobalFVGeometryCache, false);
//! We do not store the volume variables by default
SET_BOOL_PROP(ImplicitBase, EnableGlobalVolumeVariablesCache, false);
diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh
index 32cbdb2ebe..f0431f4623 100644
--- a/dumux/porousmediumflow/1p/implicit/model.hh
+++ b/dumux/porousmediumflow/1p/implicit/model.hh
@@ -93,18 +93,21 @@ public:
int eIdx = this->elementMapper().index(element);
(*rank)[eIdx] = this->gridView_().comm().rank();
- // make sure FVElementGeometry and the volume variables are bound to the element
- this->fvGeometries_().bindElement(element);
- this->curVolVars_().bindElement(element);
+ // get the local fv geometry
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
+ this->curVolVars_().bindElement(element, fvGeometry);
- const auto& fvGeometry = this->fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
const auto& spatialParams = this->problem_().spatialParams();
auto dofIdxGlobal = scv.dofIndex();
(*p)[dofIdxGlobal] = this->curVolVars(scv).pressure();
}
+
+ // velocity output
+ //velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
}
writer.attachDofData(*p, "p", isBox);
diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh
index f463222bf0..068d2bc2df 100644
--- a/dumux/porousmediumflow/2p/implicit/model.hh
+++ b/dumux/porousmediumflow/2p/implicit/model.hh
@@ -158,11 +158,12 @@ public:
(*rank)[eIdx] = this->gridView_().comm().rank();
// make sure FVElementGeometry & vol vars are bound to the element
- this->fvGeometries_().bindElement(element);
- this->curVolVars_().bindElement(element);
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bind(element);
- const auto& fvGeometry = this->fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ this->curVolVars_().bindElement(element, fvGeometry);
+
+ for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
const auto& volVars = this->curVolVars(scv);
diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh
index 90bbff0c6c..1e42005229 100644
--- a/dumux/porousmediumflow/3p/implicit/model.hh
+++ b/dumux/porousmediumflow/3p/implicit/model.hh
@@ -144,11 +144,12 @@ public:
for (const auto& scv : fvGeometry.scvs())
{
// make sure FVElementGeometry & vol vars are bound to the element
- this->fvGeometries_().bindElement(element);
- this->curVolVars_().bindElement(element);
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
- const auto& fvGeometry = this->fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ this->curVolVars_().bindElement(element, fvGeometry);
+
+ for (auto&& scv : scvs(fvGeometry))
{
auto eIdx = scv.elementIndex();
auto dofIdxGlobal = scv.dofIndex();
diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
index 9f6c94b9c2..2a726cdc78 100644
--- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
@@ -49,6 +49,10 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
@@ -116,10 +120,13 @@ public:
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
- PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace)
+ PrimaryVariables computeFlux(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const BoundaryTypes& bcTypes)
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvFace);
+ fluxVars.initAndComputeFluxes(this->problem(), element, fvGeometry, scvf);
// get upwind weights into local scope
auto massWeight = massWeight_;
diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh
index 514460534e..1f6a00d18b 100644
--- a/dumux/porousmediumflow/3p3c/implicit/model.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/model.hh
@@ -160,11 +160,12 @@ public:
for (const auto& element : elements(this->problem_().gridView()))
{
// make sure FVElementGeometry & vol vars are bound to the element
- this->fvGeometries_().bindElement(element);
- this->curVolVars_().bindElement(element);
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
- const auto& fvGeometry = this->fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ this->curVolVars_().bindElement(element, fvGeometry);
+
+ for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
if (priVarSwitch_().wasSwitched(dofIdxGlobal))
@@ -302,11 +303,12 @@ public:
(*rank)[eIdx] = this->gridView_().comm().rank();
// make sure FVElementGeometry & vol vars are bound to the element
- this->fvGeometries_().bindElement(element);
- this->curVolVars_().bindElement(element);
+ auto fvGeometry = localView(this->globalFvGeometry());
+ fvGeometry.bindElement(element);
+
+ this->curVolVars_().bindElement(element, fvGeometry);
- const auto& fvGeometry = this->fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
const auto& volVars = this->curVolVars(scv);
int dofIdxGlobal = scv.dofIndex();
diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
index 2ddd9534c0..941e763097 100644
--- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
@@ -58,10 +58,10 @@ public:
for (const auto& element : elements(problem.gridView()))
{
// make sure FVElementGeometry is bound to the element
- problem.model().fvGeometries_().bindElement(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
phasePresence_[dofIdxGlobal] = problem.initialPhasePresence(scv);
@@ -124,11 +124,12 @@ public:
for (const auto& element : elements(problem.gridView()))
{
// make sure FVElementGeometry is bound to the element
- problem.model().fvGeometries_().bindElement(element);
- volVarsVector.bindElement(element);
+ auto fvGeometry = localView(problem.model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- const auto& fvGeometry = problem.model().fvGeometries(element);
- for (const auto& scv : scvs(fvGeometry))
+ volVarsVector.bindElement(element, fvGeometry);
+
+ for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
if (!visited_[dofIdxGlobal])
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index ee999d5aae..c2fdf8a451 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -38,14 +38,18 @@ class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
- typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
// first index for the mass balance
enum { conti0EqIdx = Indices::conti0EqIdx };
@@ -95,10 +99,16 @@ public:
* \brief Evaluate the mass flux over a face of a sub control volume
* \param scvf The sub control volume face to compute the flux on
*/
- PrimaryVariables computeFlux(const SubControlVolumeFace& scvf)
+ PrimaryVariables computeFlux(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf,
+ const BoundaryTypes& bcTypes)
{
FluxVariables fluxVars;
- fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvf);
+ fluxVars.initAndComputeFluxes(this->problem(),
+ element,
+ fvGeometry,
+ scvf);
// copy weight to local scope for use in lambda expression
auto w = upwindWeight_;
@@ -114,7 +124,6 @@ public:
auto eqIdx = conti0EqIdx + phaseIdx;
flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule);
}
-
return flux;
}
diff --git a/dumux/porousmediumflow/implicit/fluxvariables.hh b/dumux/porousmediumflow/implicit/fluxvariables.hh
index a528cd9872..4eb29b8ad8 100644
--- a/dumux/porousmediumflow/implicit/fluxvariables.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariables.hh
@@ -56,6 +56,7 @@ class PorousMediumFluxVariables<TypeTag, true, false, false> : public FluxVariab
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Stencil = std::vector<IndexType>;
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
@@ -69,17 +70,22 @@ public:
void initAndComputeFluxes(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace &scvFace)
{
- ParentType::init(problem, element, scvFace);
+ ParentType::init(problem, element, fvGeometry, scvFace);
}
template<typename FunctionType>
- Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction)
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
{
- Scalar flux = AdvectionType::flux(this->problem(), this->element(), this->scvFace(), phaseIdx);
+ Scalar flux = AdvectionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->scvFace(),
+ phaseIdx);
- const auto& insideScv = this->problem().model().fvGeometries().subControlVolume(this->scvFace().insideScvIdx());
+ const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
const auto& insideVolVars = this->problem().model().curVolVars(insideScv);
const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx());
@@ -89,8 +95,11 @@ public:
return flux*upwindFunction(insideVolVars, outsideVolVars);
}
- Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
- { return AdvectionType::stencil(problem, element, scvFace); }
+ Stencil computeStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ { return AdvectionType::stencil(problem, element, fvGeometry, scvFace); }
};
@@ -109,6 +118,7 @@ class PorousMediumFluxVariables<TypeTag, true, true, false> : public FluxVariabl
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
enum
{
@@ -120,19 +130,23 @@ public:
void initAndComputeFluxes(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace &scvFace)
{
- ParentType::init(problem, element, scvFace);
+ ParentType::init(problem, element, fvGeometry, scvFace);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- advectiveFluxes_[phaseIdx] = AdvectionType::flux(problem, element, scvFace, phaseIdx);
+ advectiveFluxes_[phaseIdx] = AdvectionType::flux(problem, element, fvGeometry, scvFace, phaseIdx);
}
- Stencil computeStencil(const Problem& problem, const Element& element, const SubControlVolumeFace& scvFace)
+ Stencil computeStencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
{
// unifiy advective and diffusive stencil
- Stencil stencil = AdvectionType::stencil(problem, element, scvFace);
- Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, scvFace);
+ Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace);
+ Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, element, fvGeometry, scvFace);
stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end());
std::sort(stencil.begin(), stencil.end());
@@ -142,7 +156,7 @@ public:
}
template<typename FunctionType>
- Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction)
+ Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
{
const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx());
const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx());
@@ -155,7 +169,11 @@ public:
Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx)
{
- Scalar flux = MolecularDiffusionType::flux(this->problem(), this->element(), this->scvFace(), phaseIdx, compIdx);
+ Scalar flux = MolecularDiffusionType::flux(this->problem(),
+ this->element(),
+ this->fvGeometry(),
+ this->scvFace(),
+ phaseIdx, compIdx);
return flux;
}
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index 3bf64405ee..73c307b391 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -45,6 +45,7 @@ class PorousMediumFluxVariablesCache<TypeTag, typename std::enable_if<GET_PROP_V
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using Element = typename GridView::template Codim<0>::Entity;
@@ -52,30 +53,30 @@ class PorousMediumFluxVariablesCache<TypeTag, typename std::enable_if<GET_PROP_V
using Stencil = std::vector<IndexType>;
using TransmissibilityVector = std::vector<IndexType>;
- typedef typename GridView::ctype CoordScalar;
+ using CoordScalar = typename GridView::ctype;
static const int dim = GridView::dimension;
- typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> FeCache;
- typedef typename FeCache::FiniteElementType::Traits::LocalBasisType FeLocalBasis;
- typedef typename FeLocalBasis::Traits::JacobianType ShapeJacobian;
- typedef typename Dune::FieldVector<Scalar, 1> ShapeValue;
- typedef typename Element::Geometry::JacobianInverseTransposed JacobianInverseTransposed;
+
+ using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+ using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+ using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
+ using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
+ using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
public:
// stores required data for the flux calculation
struct FaceData
{
- std::vector<ShapeJacobian> localJacobian;
- std::vector<ShapeValue> shapeValues;
+ std::vector<ShapeJacobian> shapeJacobian;
+ std::vector<ShapeValue> shapeValue;
JacobianInverseTransposed jacInvT;
};
void update(const Problem& problem,
const Element& element,
- const typename Element::Geometry& geometry,
- const FeLocalBasis& localBasis,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace &scvFace)
{
- faceData_ = AdvectionType::calculateFaceData(problem, element, geometry, localBasis, scvFace);
+ faceData_ = AdvectionType::calculateFaceData(problem, element, fvGeometry, scvFace);
// The stencil info is obsolete for the box method.
// It is here for compatibility with cc methods
@@ -103,6 +104,7 @@ class PorousMediumFluxVariablesCache<TypeTag, typename std::enable_if<GET_PROP_V
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using Element = typename GridView::template Codim<0>::Entity;
@@ -112,11 +114,12 @@ class PorousMediumFluxVariablesCache<TypeTag, typename std::enable_if<GET_PROP_V
public:
void update(const Problem& problem,
const Element& element,
+ const FVElementGeometry& fvGeometry,
const SubControlVolumeFace &scvFace)
{
FluxVariables fluxVars;
- stencil_ = fluxVars.computeStencil(problem, element, scvFace);
- tij_ = AdvectionType::calculateTransmissibilities(problem, element, scvFace);
+ stencil_ = fluxVars.computeStencil(problem, element, fvGeometry, scvFace);
+ tij_ = AdvectionType::calculateTransmissibilities(problem, element, fvGeometry, scvFace);
}
const Stencil& stencil() const
diff --git a/test/geomechanics/elastic/elasticmatrixproblem.hh b/test/geomechanics/elastic/elasticmatrixproblem.hh
index f980d96466..612599adac 100644
--- a/test/geomechanics/elastic/elasticmatrixproblem.hh
+++ b/test/geomechanics/elastic/elasticmatrixproblem.hh
@@ -189,7 +189,7 @@ public:
PrimaryVariables values(0.0);
// inside scv
- const auto& scv = this->model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
+ auto&& scv = this->model().fvGeometries().subControlVolume(scvFace.insideScvIdx());
// get Lame parameters
Scalar lambda = this->spatialParams().lameParams(element, scv)[0];
diff --git a/test/implicit/test_fvelementgeometry.cc b/test/implicit/test_fvelementgeometry.cc
index d3e35162c3..ab4e30f4e0 100644
--- a/test/implicit/test_fvelementgeometry.cc
+++ b/test/implicit/test_fvelementgeometry.cc
@@ -32,7 +32,7 @@
#include <dune/grid/common/mcmgmapper.hh>
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
-#include <dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh>
+#include <dumux/implicit/cellcentered/tpfa/globalfvgeometry.hh>
#include <dumux/implicit/fvelementgeometry.hh>
#include <dumux/implicit/subcontrolvolume.hh>
#include <dumux/implicit/subcontrolvolumeface.hh>
@@ -119,7 +119,7 @@ int main (int argc, char *argv[]) try
if(0 != testForwardIterator(range.begin(), range.end(), op))
DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
std::cout << "-- scv center at: " << scv.center() << std::endl;
}
@@ -129,7 +129,7 @@ int main (int argc, char *argv[]) try
if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
- for (const auto& scvf : scvfs(fvGeometry))
+ for (auto&& scvf : scvfs(fvGeometry))
{
std::cout << "-- scvf center at: " << scvf.center();
if (scvf.boundary()) std::cout << " (on boundary).";
diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
index c0773b0b91..92fb753e7e 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
@@ -344,10 +344,10 @@ public:
for (const auto& element : elements(this->gridView()))
{
// make sure the FVElementGeometry is bound to the element
- this->model().fvGeometries_().bindElement(element);
- const auto& fvGeometry = this->model().fvGeometries(element);
+ auto fvGeometry = localView(this->model().globalFvGeometry());
+ fvGeometry.bindElement(element);
- for (const auto& scv : scvs(fvGeometry))
+ for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
(*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(scv);
--
GitLab