diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh
index 426567da3aad66fd288a8c3c081a2d00d9acca78..f636d72f91729867286bf38f2947596da833e234 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 0000000000000000000000000000000000000000..75e879c259df597e64b19c50d52f96cc1773c0ac
--- /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 61152459b052397766e1553f714a50909a84083a..093ba8ddf22089cf01be6049b584545d82a3fb2e 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 b62bc5e92c7d9e5a3357f8d780c61f725a0e5768..65168271c8c028e959d5067752d3bfd9ae067673 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 0000000000000000000000000000000000000000..5d843b48095c63553331fb9bcc3bfc235760d2d3
--- /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 48c11d229f87c976f351d25437dc735b2ae6de00..0000000000000000000000000000000000000000
--- 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 0000000000000000000000000000000000000000..e1efe5803a13509de9a0f488bd71439f290ae88b
--- /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 afe38dc3f8b2023eeef0069bbec5c06d3d2618fe..489dead98995b32bdfed9bcdd1c301b77874896c 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 df886d74d351a46a4a3e58b029edd245d67efbda..410a2f272a57d5638a107be1e03ae09324df7c26 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 b96d5085c58f0d20dbcac3e393b966cbb64aeb77..1f297a77e203c8aed3fb21421b304540454dfeb9 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 313f39f796eee3fc2b9541c676baf366f2def39b..31e88237bf17c0e26d5841a4ad0588e3f2e8a48b 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 859f16d12f3ceea8bbf5262fb210f4891ac7c24c..a44c8084f4a6b04cbf442ed3aa3875f183721656 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 98c0a346d8b9048508884a9963131ddf5cacec57..ef2b34e86caea8e681e1d48b82563c982d3f0fbc 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 ff7ebeb2ba8ab4c57db241f4588d5f05ef4d7663..9bb40dd86338a59dfac0202c5efd7bbc3b3e9443 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 c24d224618ad2a754c52a16ac3c651c5ae1d71c0..345e94d58c6ca6c28d08201b567c52643ea1bfee 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 cddca62731be7c0ca1fb190a13b38dd7a785f4ed..6f6d5cc37b4debe234faa92ba53c4516d2e41056 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 0000000000000000000000000000000000000000..9acef62f15faa0e654d11da4237d68b541a964d9
--- /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 1a239b4e38eee895288e0a3940b2ba40aac6250f..35e3e1f15dd1056cf95757105513a4287dd55455 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 aaba6ee555717cba8bf670ebffc7a45bb9bc9222..40d78cfaee6a1fba9e7b0d3d9d338bfbdc28fc19 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 9a3d763c5c3d3cca5de94071ee97353608e2d30d..6c43bc501528d90359867f85857cc17ba500dba4 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 39b62356521f3c019af01fc8675799e21a315880..1fcc38cf146981f41c2e35a7a3a681ecdecea2d7 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 de8061254d48ce12d473e717c04a6ee5430b64b1..52eb6dfd7204cf0337a5cbb09d605ba8e744884f 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 7640c0be3f4d01570134921b1e5649be370fcb9b..0000000000000000000000000000000000000000
--- 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 0000000000000000000000000000000000000000..81c01a1ab86a613a978a8e83fff03b821a322b00
--- /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 56f141514ad036325286abb7708a21c676888f62..15722ac3e3912547cad7071f3fe1f84b7c80740e 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 ab732275d7e7c33d073146635603ee9b5d627879..01e962c3153a477a82e017823449532c9869603f 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 809b59ffc4de898c014384c13b2443791a01eafa..b2c5c13b1bbe0cccaabc0ddc28e61919dfcb6152 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 09e16ae56f1ea0a3541e104ac21b831cf252ef5f..e824559e463b2648791675fba78c8fffee1a2601 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 69ec63c1c01fbee0a66ef3db3cf9b689c4f2e15d..1a8e24d07213544a52c16c9ebd89c98cc9b65db9 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 bd1b8af5bcab8ef30229aae3959d50926903642f..38ef7f2044779171b49a98b01f9741cf454ba947 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 0f20ee7e1c45e71d9f0da2ef414ec1c91b9e22e7..d8ff00658203fa5c2a11a6e5f0169f4006484578 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 3681ef036c5559759336bcbc84057b2c62819ed6..71d919fcea438ad0debbdb6e8e070dd8200b2b8f 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 a11630ed96275895801ba84daf099c7f445402c6..addff4924c4379d128d9294ae79b5a4ace7f7e60 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 a38fe2d64063fd7eb54d27509a5209cdf0805941..3ddce0a1cbefeb8e6d1e263e6933c629cafb9603 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 4766debb5ff239e8b5de58517485275aac1c8627..2882022e1db9af03bb91e7869f5816a8538cddca 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 cc904e2b56783dc6a1d67589904e457e8811a284..5ce873a384f1fe92c5ec7c18afba04016c6416df 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 96c8ee5cc3df1a31bc9c2def71f3be70cede90e0..2396859acdb509a28008976e8220800f28f70e3e 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 8a76072f7185cf7a1e1350ef066b23f792a3fcb7..4088bf1cc96857682707a7adcaadb523044012ab 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 77c9cd298748e878be0b4a5c12cf4de12bcfa66d..916cb385dd3cdf2fea61a9108cb78679a5a61bb2 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 c4effc39f902f57720a44bba9ea537cd1dd40817..cff3fdf63c97ca79a8bcde659272416b4fb40e2c 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 2d4219e401879e44f97545bcfb3f40e5e81f5644..864a2b02027ed141261ac6b0a6b5534ec64411e2 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 32cbdb2ebe75673526392457b0b631d65888a8e1..f0431f462395e00304d0dbb2a7e490b70ffbd436 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 f463222bf0cad0922a09b0e3a919ea0c524a3d2f..068d2bc2df15f8f12efac6a5b1f9d5742b8a0ede 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 90bbff0c6c86d8f9ed0312317efb4b5b367f70e9..1e42005229172d25294a326fd8bb5cbca8ad94f6 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 9f6c94b9c208d590cf3497fd0379a8386d965e47..2a726cdc787503bc8c40db21ba7fe637d9249664 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 514460534e83557f5f4a4b221ba64dc0f065af3f..1f6a00d18b568ceb971dfdddde8de0c7d0bcdb8a 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 2ddd9534c0b89e50e1358bfc910ea1d60bc6b9b8..941e76309767c242df990ccb20911af703c60f93 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 ee999d5aaedbaf090add9ff343cde9738be2a992..c2fdf8a451fa7691c82fcf57d9891ee97536c7d1 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 a528cd9872120c7b479e97c192de8befafdb2fbd..4eb29b8ad8f7f237152c150fd037ff296cb8d96a 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 3bf64405eed1b750f9d905c1a3a7d5b22146db50..73c307b3915d86a3e6004320bb9fc969128d90db 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 f980d9646670e541b23e682f830207e74d07c7ff..612599adacde51b4c874ba96f5e7101975854405 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 d3e35162c3e4532378b9e3386a819132d2750aff..ab4e30f4e0367c66b979db651a5c104b6fc1accb 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 c0773b0b9146de9b2860c918a158bfd11fae851c..92fb753e7ed91ba28f68be7bd365ad49beafe081 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);