From b7fbfc28023ad1f3062de0a9405691469677656b Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Fri, 11 Nov 2016 15:35:22 +0100
Subject: [PATCH] [staggeredGrid] Add basic folder structure and first test
---
dumux/discretization/staggered/darcyslaw.hh | 226 +++++++++++
.../staggered/fvelementgeometry.hh | 355 +++++++++++++++++
.../staggered/globalfluxvariablescache.hh | 138 +++++++
.../staggered/globalfvgeometry.hh | 363 ++++++++++++++++++
.../staggered/staggeredgeometryhelper.hh | 39 ++
.../staggered/subcontrolvolume.hh | 139 +++++++
.../staggered/subcontrolvolumeface.hh | 220 +++++++++++
dumux/implicit/staggered/properties.hh | 52 +++
dumux/implicit/staggered/propertydefaults.hh | 115 ++++++
test/discretization/CMakeLists.txt | 1 +
test/discretization/staggered/CMakeLists.txt | 2 +
.../staggered/test_staggered.cc | 149 +++++++
12 files changed, 1799 insertions(+)
create mode 100644 dumux/discretization/staggered/darcyslaw.hh
create mode 100644 dumux/discretization/staggered/fvelementgeometry.hh
create mode 100644 dumux/discretization/staggered/globalfluxvariablescache.hh
create mode 100644 dumux/discretization/staggered/globalfvgeometry.hh
create mode 100644 dumux/discretization/staggered/staggeredgeometryhelper.hh
create mode 100644 dumux/discretization/staggered/subcontrolvolume.hh
create mode 100644 dumux/discretization/staggered/subcontrolvolumeface.hh
create mode 100644 dumux/implicit/staggered/properties.hh
create mode 100644 dumux/implicit/staggered/propertydefaults.hh
create mode 100644 test/discretization/staggered/CMakeLists.txt
create mode 100644 test/discretization/staggered/test_staggered.cc
diff --git a/dumux/discretization/staggered/darcyslaw.hh b/dumux/discretization/staggered/darcyslaw.hh
new file mode 100644
index 0000000000..06c5d58d1c
--- /dev/null
+++ b/dumux/discretization/staggered/darcyslaw.hh
@@ -0,0 +1,226 @@
+// -*- 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 This file contains the data which is required to calculate
+ * volume and mass fluxes of fluid phases over a face of a finite volume by means
+ * of the Darcy approximation. Specializations are provided for the different discretization methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_DARCYS_LAW_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_DARCYS_LAW_HH
+
+#include <memory>
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/methods.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+// forward declaration of properties
+NEW_PROP_TAG(ProblemEnableGravity);
+}
+
+/*!
+ * \ingroup DarcysLaw
+ * \brief Specialization of Darcy's Law for the CCTpfa method. TODO: dummy, remove this class
+ */
+template <class TypeTag>
+class DarcysLawImplementation<TypeTag, DiscretizationMethods::Staggered>
+{
+ 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);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FluxVarsCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+
+ 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} ;
+
+ 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 ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvFace,
+ int phaseIdx,
+ const FluxVarsCache& fluxVarsCache)
+ {
+ const auto& tij = fluxVarsCache.tij();
+
+ // Get the inside and outside volume variables
+ const auto& insideScv = fvGeometry.scv(scvFace.insideScvIdx());
+ const auto& insideVolVars = elemVolVars[insideScv];
+ const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
+
+ auto hInside = insideVolVars.pressure(phaseIdx);
+ auto hOutside = outsideVolVars.pressure(phaseIdx);
+
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+ {
+ // do averaging for the density
+ const auto rhoInside = insideVolVars.density(phaseIdx);
+ const auto rhoOutide = outsideVolVars.density(phaseIdx);
+ const auto rho = (rhoInside + rhoOutide)*0.5;
+
+ // ask for the gravitational acceleration in the inside neighbor
+ const auto xInside = insideScv.center();
+ const auto gInside = problem.gravityAtPos(xInside);
+
+ hInside -= rho*(gInside*xInside);
+
+ // and the outside neighbor
+ if (scvFace.boundary())
+ {
+ const auto xOutside = scvFace.center();
+ const auto gOutside = problem.gravityAtPos(xOutside);
+ hOutside -= rho*(gOutside*xOutside);
+ }
+ else
+ {
+ const auto outsideScvIdx = scvFace.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);
+ }
+ }
+
+ return tij*(hInside - hOutside);
+ }
+
+ static Stencil stencil(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvFace)
+ {
+ Stencil stencil;
+ if (!scvFace.boundary())
+ {
+ stencil.push_back(scvFace.insideScvIdx());
+ stencil.push_back(scvFace.outsideScvIdx());
+ }
+ else
+ stencil.push_back(scvFace.insideScvIdx());
+
+ return stencil;
+ }
+
+ // 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 FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvFace)
+ {
+ Scalar tij;
+
+ const auto insideScvIdx = scvFace.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+ const auto& insideVolVars = elemVolVars[insideScvIdx];
+ const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv, insideVolVars);
+ Scalar ti = calculateOmega_(problem, scvFace, insideK, element, insideScv);
+
+ if (!scvFace.boundary())
+ {
+ const auto outsideScvIdx = scvFace.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& outsideVolVars = elemVolVars[outsideScvIdx];
+ const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv, outsideVolVars);
+ Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideElement, outsideScv);
+
+ // check for division by zero!
+ if (ti*tj <= 0.0)
+ tij = 0;
+ else
+ tij = scvFace.area()*(ti * tj)/(ti + tj);
+ }
+ else
+ {
+ tij = scvFace.area()*ti;
+ }
+
+ return tij;
+ }
+
+private:
+
+ 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);
+
+ auto distanceVector = scvFace.center();
+ distanceVector -= scv.center();
+ distanceVector /= distanceVector.two_norm2();
+
+ Scalar omega = Knormal * distanceVector;
+ omega *= problem.boxExtrusionFactor(element, scv);
+
+ return omega;
+ }
+
+ 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();
+ distanceVector /= distanceVector.two_norm2();
+
+ Scalar omega = K * (distanceVector * scvFace.unitOuterNormal());
+ omega *= problem.boxExtrusionFactor(element, scv);
+
+ return omega;
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
new file mode 100644
index 0000000000..7e12e9fc43
--- /dev/null
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -0,0 +1,355 @@
+// -*- 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 staggered 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_STAGGERED_FV_ELEMENT_GEOMETRY_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_FV_ELEMENT_GEOMETRY_HH
+
+#include <dune/common/iteratorrange.hh>
+
+#include <dumux/discretization/scvandscvfiterators.hh>
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/discretization/staggered/subcontrolvolume.hh>
+
+namespace Dumux
+{
+
+//! forward declaration of the global finite volume geometry
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class StaggeredGlobalFVGeometry;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for staggered models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class StaggeredFVElementGeometry
+{};
+
+//! specialization in case the FVElementGeometries are stored globally
+//! In this case we just forward internally to the global object
+template<class TypeTag>
+class StaggeredFVElementGeometry<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 = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
+ using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
+
+public:
+ //! Constructor
+ StaggeredFVElementGeometry(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 StaggeredFVElementGeometry& 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 StaggeredFVElementGeometry& 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 StaggeredFVElementGeometry<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 = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
+ using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
+
+public:
+ //! Constructor
+ StaggeredFVElementGeometry(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))
+ {
+ neighborScvs_.reserve(element.subEntities(1));
+ neighborScvfIndices_.reserve(element.subEntities(1));
+ 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,
+ intersection.geometry(),
+ scvFaceIndices[scvfCounter],
+ std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
+ );
+
+ scvfIndices_.emplace_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,
+ intersection.geometry(),
+ scvFaceIndices[scvfCounter],
+ std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]})
+ );
+
+ neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
+ }
+ // always increase local counter
+ scvfCounter++;
+ }
+ else if (intersection.boundary())
+ {
+ 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/staggered/globalfluxvariablescache.hh b/dumux/discretization/staggered/globalfluxvariablescache.hh
new file mode 100644
index 0000000000..8cd7bc756a
--- /dev/null
+++ b/dumux/discretization/staggered/globalfluxvariablescache.hh
@@ -0,0 +1,138 @@
+// -*- 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 The global object of flux var caches
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_GLOBAL_FLUXVARSCACHE_HH
+#define DUMUX_DISCRETIZATION_CC_GLOBAL_FLUXVARSCACHE_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/cellcentered/elementfluxvariablescache.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the flux variables cache vector, we store one cache per face
+ */
+template<class TypeTag, bool EnableGlobalFluxVariablesCache>
+class CCGlobalFluxVariablesCache;
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Spezialization when caching globally
+ */
+template<class TypeTag>
+class CCGlobalFluxVariablesCache<TypeTag, true>
+{
+ // the local class needs access to the problem
+ friend CCElementFluxVariablesCache<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 FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+public:
+ // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
+ void update(Problem& problem)
+ {
+ problemPtr_ = &problem;
+ 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
+ auto fvGeometry = localView(globalFvGeometry);
+ fvGeometry.bind(element);
+
+ auto elemVolVars = localView(problem.model().curGlobalVolVars());
+ elemVolVars.bind(element, fvGeometry, problem.model().curSol());
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf);
+ }
+ }
+ }
+
+ /*!
+ * \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 inline ElementFluxVariablesCache localView(const CCGlobalFluxVariablesCache& global)
+ { return ElementFluxVariablesCache(global); }
+
+private:
+ // access operators in the case of caching
+ const FluxVariablesCache& operator [](IndexType scvfIdx) const
+ { return fluxVarsCache_[scvfIdx]; }
+
+ FluxVariablesCache& operator [](IndexType scvfIdx)
+ { return fluxVarsCache_[scvfIdx]; }
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+
+ std::vector<FluxVariablesCache> fluxVarsCache_;
+ std::vector<IndexType> globalScvfIndices_;
+};
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Spezialization when not using global caching
+ */
+template<class TypeTag>
+class CCGlobalFluxVariablesCache<TypeTag, false>
+{
+ // the local class needs access to the problem
+ friend CCElementFluxVariablesCache<TypeTag, false>;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+
+public:
+ // When global flux variables caching is disabled, we don't need to update the cache
+ void update(Problem& problem)
+ { problemPtr_ = &problem; }
+
+ /*!
+ * \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 inline ElementFluxVariablesCache localView(const CCGlobalFluxVariablesCache& global)
+ { return ElementFluxVariablesCache(global); }
+
+private:
+
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/globalfvgeometry.hh b/dumux/discretization/staggered/globalfvgeometry.hh
new file mode 100644
index 0000000000..bcfb8d03c1
--- /dev/null
+++ b/dumux/discretization/staggered/globalfvgeometry.hh
@@ -0,0 +1,363 @@
+// -*- 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 staggered models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element of the grid partition.
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FVGEOMETRY_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_FVGEOMETRY_HH
+
+#include <dumux/common/elementmap.hh>
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/discretization/staggered/fvelementgeometry.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the finite volume geometry vector for staggered models
+ * This builds up the sub control volumes and sub control volume faces
+ * for each element.
+ */
+template<class TypeTag, bool EnableGlobalFVGeometryCache>
+class StaggeredGlobalFVGeometry
+{};
+
+// specialization in case the FVElementGeometries are stored globally
+template<class TypeTag>
+class StaggeredGlobalFVGeometry<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 and the fv element geometry
+ //! as they are globally cached
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+public:
+ //! Constructor
+ StaggeredGlobalFVGeometry(const GridView& gridView)
+ : gridView_(gridView), elementMap_(gridView) {}
+
+ //! 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();
+ }
+
+ //! The total number of boundary sub control volume faces
+ 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.elementIndex()); }
+
+ // 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;
+
+ // clear containers (necessary after grid refinement)
+ scvs_.clear();
+ scvfs_.clear();
+ scvfIndicesOfScv_.clear();
+ elementMap_.clear();
+
+ // determine size of containers
+ const int numElements = gridView_.size(0);
+ IndexType numScvs = gridView_.size(0);
+ IndexType numScvf = 0;
+ for (const auto& element : elements(gridView_))
+ numScvf += element.subEntities(1);
+
+ // reserve memory
+ elementMap_.resize(numScvs);
+ scvs_.resize(numScvs);
+ scvfs_.reserve(numScvf);
+ scvfIndicesOfScv_.resize(numScvs);
+
+ // Build the scvs and scv faces
+ IndexType scvfIdx = 0;
+ numBoundaryScvf_ = 0;
+ for (const auto& element : elements(gridView_))
+ {
+ auto eIdx = problem.elementMapper().index(element);
+ scvs_[eIdx] = SubControlVolume(element.geometry(), eIdx);
+
+ // fill the element map with seeds
+ elementMap_[eIdx] = element.seed();
+
+ // the element-wise index sets for finite volume geometry
+ std::vector<IndexType> scvfsIndexSet;
+ scvfsIndexSet.reserve(element.subEntities(1));
+ for (const auto& intersection : intersections(gridView_, element))
+ {
+ //TODO: use proper intersection mapper!
+ const auto inIdx = intersection.indexInInside();
+ const auto globalIsIdx = gridView_.indexSet().subIndex(element, inIdx, dim-1) + numElements;
+
+ int oppoSiteIdx;
+
+ if(inIdx % 2) // face index is odd
+ {
+ oppoSiteIdx = inIdx -1;
+ }
+ else
+ {
+ oppoSiteIdx = inIdx +1;
+ }
+
+// std::cout << "faceSelf: " << inIdx << ", oppo: " << oppoSiteIdx << std::endl;
+ // inner sub control volume faces
+ if (intersection.neighbor())
+ {
+ auto nIdx = problem.elementMapper().index(intersection.outside());
+ scvfs_.emplace_back(intersection,
+ intersection.geometry(),
+ scvfIdx,
+ std::vector<IndexType>({eIdx, nIdx}),
+ globalIsIdx
+ );
+ scvfsIndexSet.push_back(scvfIdx++);
+ }
+ // boundary sub control volume faces
+ else if (intersection.boundary())
+ {
+ scvfs_.emplace_back(intersection,
+ intersection.geometry(),
+ scvfIdx,
+ std::vector<IndexType>({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
+ globalIsIdx
+ );
+ scvfsIndexSet.push_back(scvfIdx++);
+ }
+ }
+
+ // Save the scvf indices belonging to this scv to build up fv element geometries fast
+ scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
+ }
+ }
+
+ /*!
+ * \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 inline FVElementGeometry localView(const StaggeredGlobalFVGeometry& 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];
+ }
+
+private:
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+
+ const GridView gridView_;
+ Dumux::ElementMap<GridView> elementMap_;
+ std::vector<SubControlVolume> scvs_;
+ std::vector<SubControlVolumeFace> scvfs_;
+ std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
+ IndexType numBoundaryScvf_;
+};
+
+// specialization in case the FVElementGeometries are not stored
+template<class TypeTag>
+class StaggeredGlobalFVGeometry<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 Element = typename GridView::template Codim<0>::Entity;
+ //! The local fvGeometry needs access to the problem
+ friend typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+
+public:
+ //! Constructor
+ StaggeredGlobalFVGeometry(const GridView& gridView)
+ : gridView_(gridView), elementMap_(gridView) {}
+
+ //! The total number of sub control volumes
+ std::size_t numScv() const
+ {
+ return numScvs_;
+ }
+
+ //! The total number of sub control volume faces
+ std::size_t numScvf() const
+ {
+ return numScvf_;
+ }
+
+ //! The total number of boundary sub control volume faces
+ 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.elementIndex()); }
+
+ // 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();
+
+ // reserve memory or resize the containers
+ numScvs_ = gridView_.size(0);
+ numScvf_ = 0;
+ numBoundaryScvf_ = 0;
+ elementMap_.resize(numScvs_);
+ scvfIndicesOfScv_.resize(numScvs_);
+ neighborVolVarIndices_.resize(numScvs_);
+
+ // Build the SCV and SCV face
+ for (const auto& element : elements(gridView_))
+ {
+ auto eIdx = problem.elementMapper().index(element);
+
+ // fill the element map with seeds
+ elementMap_[eIdx] = element.seed();
+
+ // the element-wise index sets for finite volume geometry
+ auto numLocalFaces = element.subEntities(1);
+ std::vector<IndexType> scvfsIndexSet;
+ std::vector<IndexType> neighborVolVarIndexSet;
+ scvfsIndexSet.reserve(numLocalFaces);
+ neighborVolVarIndexSet.reserve(numLocalFaces);
+ for (const auto& intersection : intersections(gridView_, element))
+ {
+ // inner sub control volume faces
+ if (intersection.neighbor())
+ {
+ scvfsIndexSet.push_back(numScvf_++);
+ neighborVolVarIndexSet.push_back(problem.elementMapper().index(intersection.outside()));
+ }
+ // boundary sub control volume faces
+ else if (intersection.boundary())
+ {
+ scvfsIndexSet.push_back(numScvf_++);
+ neighborVolVarIndexSet.push_back(numScvs_ + numBoundaryScvf_++);
+ }
+ }
+
+ // store the sets of indices in the data container
+ scvfIndicesOfScv_[eIdx] = scvfsIndexSet;
+ neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
+ }
+ }
+
+ const std::vector<IndexType>& scvfIndicesOfScv(IndexType scvIdx) const
+ { return scvfIndicesOfScv_[scvIdx]; }
+
+ const std::vector<IndexType>& neighborVolVarIndices(IndexType scvIdx) const
+ { return neighborVolVarIndices_[scvIdx]; }
+
+ /*!
+ * \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 inline FVElementGeometry localView(const StaggeredGlobalFVGeometry& global)
+ { return FVElementGeometry(global); }
+
+private:
+ const Problem& problem_() const
+ { return *problemPtr_; }
+
+ const Problem* problemPtr_;
+
+ const GridView gridView_;
+
+ // Information on the global number of geometries
+ IndexType numScvs_;
+ IndexType numScvf_;
+ IndexType numBoundaryScvf_;
+
+ // vectors that store the global data
+ Dumux::ElementMap<GridView> elementMap_;
+ std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
+ std::vector<std::vector<IndexType>> neighborVolVarIndices_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/staggeredgeometryhelper.hh b/dumux/discretization/staggered/staggeredgeometryhelper.hh
new file mode 100644
index 0000000000..6d7339bbde
--- /dev/null
+++ b/dumux/discretization/staggered/staggeredgeometryhelper.hh
@@ -0,0 +1,39 @@
+// -*- 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 staggered discretizazion method
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/staggered/subcontrolvolume.hh b/dumux/discretization/staggered/subcontrolvolume.hh
new file mode 100644
index 0000000000..6cce08266d
--- /dev/null
+++ b/dumux/discretization/staggered/subcontrolvolume.hh
@@ -0,0 +1,139 @@
+// -*- 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 sub control volume
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_SUBCONTROLVOLUME_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_SUBCONTROLVOLUME_HH
+
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumebase.hh>
+#include <dumux/common/optional.hh>
+
+namespace Dumux
+{
+template<class G, typename I>
+class StaggeredSubControlVolume : public SubControlVolumeBase<StaggeredSubControlVolume<G, I>, G, I>
+{
+ using ParentType = SubControlVolumeBase<StaggeredSubControlVolume<G, I>, G, I>;
+ using Geometry = G;
+ using IndexType = I;
+
+ using Scalar = typename Geometry::ctype;
+ enum { dimworld = Geometry::coorddimension };
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+public:
+ // the default constructor
+ StaggeredSubControlVolume() = default;
+
+ // the contructor in the cc case
+ StaggeredSubControlVolume(Geometry&& geometry,
+ IndexType elementIndex)
+ : ParentType(), geometry_(std::move(geometry)), elementIndex_(elementIndex) {}
+
+ //! The copy constrcutor
+ StaggeredSubControlVolume(const StaggeredSubControlVolume& other) = default;
+
+ //! The move constrcutor
+ StaggeredSubControlVolume(StaggeredSubControlVolume&& other) = default;
+
+ //! The copy assignment operator
+ StaggeredSubControlVolume& operator=(const StaggeredSubControlVolume& other)
+ {
+ // We want to use the default copy/move assignment.
+ // But since geometry is not copy assignable :( we
+ // have to construct it again
+ geometry_.release();
+ geometry_.emplace(other.geometry_.value());
+ elementIndex_ = other.elementIndex_;
+ return *this;
+ }
+
+ //! The move assignment operator
+ StaggeredSubControlVolume& operator=(StaggeredSubControlVolume&& other)
+ {
+ // We want to use the default copy/move assignment.
+ // But since geometry is not copy assignable :( we
+ // have to construct it again
+ geometry_.release();
+ geometry_.emplace(other.geometry_.value());
+ elementIndex_ = std::move(other.elementIndex_);
+ return *this;
+ }
+
+ //! The center of the sub control volume
+ GlobalPosition center() const
+ {
+ return geometry().center();
+ }
+
+ //! The volume of the sub control volume
+ Scalar volume() const
+ {
+ return geometry().volume();
+ }
+
+ //! The geometry of the sub control volume
+ // e.g. for integration
+ const Geometry& geometry() const
+ {
+ return geometry_.value();
+ }
+
+ //! The global index of this scv
+ IndexType index() const
+ {
+ return elementIndex();
+ }
+
+ //! The index of the dof this scv is embedded in
+ IndexType dofIndex() const
+ {
+ return elementIndex();
+ }
+
+ // The position of the dof this scv is embedded in
+ GlobalPosition dofPosition() const
+ {
+ return geometry().center();
+ }
+
+ //! The global index of the element this scv is embedded in
+ IndexType elementIndex() const
+ {
+ return elementIndex_;
+ }
+
+ //! Return the corner for the given local index
+ GlobalPosition corner(unsigned int localIdx) const
+ {
+ assert(localIdx < geometry().corners().size() && "provided index exceeds the number of corners");
+ return geometry().corners(localIdx);
+ }
+
+private:
+ // Work around the fact that geometry is not default constructible
+ Optional<Geometry> geometry_;
+ IndexType elementIndex_;
+};
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/subcontrolvolumeface.hh b/dumux/discretization/staggered/subcontrolvolumeface.hh
new file mode 100644
index 0000000000..a137e17722
--- /dev/null
+++ b/dumux/discretization/staggered/subcontrolvolumeface.hh
@@ -0,0 +1,220 @@
+// -*- 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 sub control volume face
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_SUBCONTROLVOLUMEFACE_HH
+
+#include <utility>
+#include <dune/common/fvector.hh>
+#include <dumux/discretization/subcontrolvolumefacebase.hh>
+#include <dumux/common/optional.hh>
+
+namespace Dumux
+{
+
+template<class G, typename I>
+class StaggeredSubFace
+{
+ using Geometry = G;
+ using IndexType = I;
+ using Scalar = typename Geometry::ctype;
+
+private:
+ std::vector<std::pair<int,int>> velocityDofIdxPair_;
+ std::vector<Scalar> distance_;
+ std::pair<int,int> elementPair_;
+ int commonVertexIdx_;
+};
+
+
+/*!
+ * \ingroup Discretization
+ * \brief Class for a sub control volume face in the box method, i.e a part of the boundary
+ * of a sub control volume we compute fluxes on. We simply use the base class here.
+ */
+template<class G, typename I>
+class StaggeredSubControlVolumeFace : public SubControlVolumeFaceBase<StaggeredSubControlVolumeFace<G, I>, G, I>
+{
+ using ParentType = SubControlVolumeFaceBase<StaggeredSubControlVolumeFace<G, I>, G, I>;
+ using Geometry = G;
+ using IndexType = I;
+
+ using Scalar = typename Geometry::ctype;
+ static const int dim = Geometry::mydimension;
+ static const int dimworld = Geometry::coorddimension;
+
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+ using LocalPosition = Dune::FieldVector<Scalar, dim>;
+
+ using StaggeredSubFace = Dumux::StaggeredSubFace<G,I>;
+
+public:
+ // the default constructor
+ StaggeredSubControlVolumeFace() = default;
+
+ //! Constructor with intersection
+ template <class Intersection>
+ StaggeredSubControlVolumeFace(const Intersection& is,
+ const typename Intersection::Geometry& isGeometry,
+ IndexType scvfIndex,
+ const std::vector<IndexType>& scvIndices,
+ const int selfIdx
+ )
+ : ParentType(),
+ geomType_(isGeometry.type()),
+ area_(isGeometry.volume()),
+ center_(isGeometry.center()),
+ unitOuterNormal_(is.centerUnitOuterNormal()),
+ scvfIndex_(scvfIndex),
+ scvIndices_(scvIndices),
+ boundary_(is.boundary()),
+ selfIdx_(selfIdx)
+ {
+ corners_.resize(isGeometry.corners());
+ for (int i = 0; i < isGeometry.corners(); ++i)
+ corners_[i] = isGeometry.corner(i);
+
+// subfaces_.resize(2);
+
+// const auto& element = is.inside();
+// const int inIdx = is.indexInInside();
+// gridView.indexSet().subIndex(element, inIdx, dimWorld-1);
+
+ }
+
+ /*//! The copy constrcutor
+ StaggeredSubControlVolumeFace(const StaggeredSubControlVolumeFace& other) = delete;
+
+ //! The move constrcutor
+ StaggeredSubControlVolumeFace(StaggeredSubControlVolumeFace&& other) = default;
+
+ //! The copy assignment operator
+ StaggeredSubControlVolumeFace& operator=(const StaggeredSubControlVolumeFace& other) = delete;
+
+ //! The move assignment operator
+ // StaggeredSubControlVolumeFace& operator=(StaggeredSubControlVolumeFace&& other)
+ // {
+ // // We want to use the default copy/move assignment.
+ // // But since geometry is not copy assignable :( we
+ // // have to construct it again
+ // geometry_.release();
+ // geometry_.emplace(other.geometry_.value());
+ // unitOuterNormal_ = std::move(other.unitOuterNormal_);
+ // scvfIndex_ = std::move(other.scvfIndex_);
+ // scvIndices_ = std::move(other.scvIndices_);
+ // boundary_ = std::move(other.boundary_);
+ // return *this;
+ // }*/
+
+ //! The center of the sub control volume face
+ GlobalPosition center() const
+ {
+ return center_;
+ }
+
+ //! The integration point for flux evaluations in global coordinates
+ GlobalPosition ipGlobal() const
+ {
+ // Return center for now
+ return center_;
+ }
+
+ //! The area of the sub control volume face
+ Scalar area() const
+ {
+ return area_;
+ }
+
+ //! returns bolean if the sub control volume face is on the boundary
+ bool boundary() const
+ {
+ return boundary_;
+ }
+
+ GlobalPosition unitOuterNormal() const
+ {
+ return unitOuterNormal_;
+ }
+
+ //! index of the inside sub control volume for spatial param evaluation
+ IndexType insideScvIdx() const
+ {
+ return scvIndices_[0];
+ }
+
+ //! index of the outside sub control volume for spatial param evaluation
+ // This results in undefined behaviour if boundary is true
+ IndexType outsideScvIdx() const
+ {
+ return scvIndices_[1];
+ }
+
+ //! The global index of this sub control volume face
+ IndexType index() const
+ {
+ return scvfIndex_;
+ }
+
+ GlobalPosition corner(unsigned int localIdx) const
+ {
+ assert(localIdx < corners_.size() && "provided index exceeds the number of corners");
+ return corners_[localIdx];
+ }
+
+ //! The geometry of the sub control volume face
+ const Geometry geometry() const
+ {
+ return Geometry(geomType_, corners_);
+ }
+
+ //! The global index of this sub control volume face
+ IndexType dofIndexSelf() const
+ {
+ return selfIdx_;
+ }
+
+ //! The global index of this sub control volume face
+ IndexType dofIndexOpposite() const
+ {
+ return oppositeIdx_;
+ }
+
+private:
+ Dune::GeometryType geomType_;
+ std::vector<GlobalPosition> corners_;
+ Scalar area_;
+ GlobalPosition center_;
+ GlobalPosition unitOuterNormal_;
+ IndexType scvfIndex_;
+ std::vector<IndexType> scvIndices_;
+ bool boundary_;
+
+ int selfIdx_;
+ int oppositeIdx_;
+ std::vector<StaggeredSubFace> subfaces_;
+};
+
+
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/staggered/properties.hh b/dumux/implicit/staggered/properties.hh
new file mode 100644
index 0000000000..73510d65ab
--- /dev/null
+++ b/dumux/implicit/staggered/properties.hh
@@ -0,0 +1,52 @@
+// -*- 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_STAGGERED_PROPERTIES_HH
+#define DUMUX_STAGGERED_PROPERTIES_HH
+
+#include <dumux/implicit/properties.hh>
+
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup StaggeredModel
+ * \file
+ * \brief Specify the shape functions, operator assemblers, etc
+ * used for the StaggeredModel.
+ */
+namespace Dumux
+{
+
+namespace Properties
+{
+// \{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for models based on staggered scheme
+NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(ImplicitBase));
+}
+}
+
+// \}
+
+#include <dumux/implicit/staggered/propertydefaults.hh>
+
+#endif
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
new file mode 100644
index 0000000000..6f49cbfaba
--- /dev/null
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -0,0 +1,115 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup CCTpfaProperties
+ * \ingroup StaggeredModel
+ * \file
+ *
+ * \brief Default properties for cell centered models
+ */
+#ifndef DUMUX_STAGGERED_PROPERTY_DEFAULTS_HH
+#define DUMUX_STAGGERED_PROPERTY_DEFAULTS_HH
+
+#include <dumux/implicit/propertydefaults.hh>
+#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
+#include <dumux/discretization/staggered/globalfvgeometry.hh>
+#include <dumux/discretization/staggered/fvelementgeometry.hh>
+#include <dumux/discretization/staggered/subcontrolvolumeface.hh>
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/stencils.hh>
+
+namespace Dumux {
+
+// forward declarations
+template<class TypeTag> class CCElementBoundaryTypes;
+
+namespace Properties {
+//! Set the corresponding discretization method property
+SET_PROP(StaggeredModel, DiscretizationMethod)
+{
+ static const DiscretizationMethods value = DiscretizationMethods::Staggered;
+};
+
+//! Set the default for the global finite volume geometry
+SET_TYPE_PROP(StaggeredModel, GlobalFVGeometry, StaggeredGlobalFVGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
+
+//! Set the default for the local finite volume geometry
+SET_TYPE_PROP(StaggeredModel, FVElementGeometry, StaggeredFVElementGeometry<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFVGeometryCache)>);
+
+//! The sub control volume
+SET_PROP(StaggeredModel, SubControlVolume)
+{
+private:
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using ScvGeometry = typename Grid::template Codim<0>::Geometry;
+ using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
+public:
+ typedef Dumux::StaggeredSubControlVolume<ScvGeometry, IndexType> type;
+};
+
+SET_PROP(StaggeredModel, SubControlVolumeFace)
+{
+private:
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using ScvfGeometry = typename Grid::template Codim<1>::Geometry;
+ using IndexType = typename Grid::LeafGridView::IndexSet::IndexType;
+public:
+ typedef Dumux::StaggeredSubControlVolumeFace<ScvfGeometry, IndexType> type;
+};
+
+//! Set the default for the ElementBoundaryTypes
+SET_TYPE_PROP(StaggeredModel, ElementBoundaryTypes, Dumux::CCElementBoundaryTypes<TypeTag>);
+
+//! Mapper for the degrees of freedoms.
+SET_TYPE_PROP(StaggeredModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper));
+
+//! The global current volume variables vector class
+SET_TYPE_PROP(StaggeredModel, GlobalVolumeVariables, Dumux::CCGlobalVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
+
+//! The global flux variables cache vector class
+SET_TYPE_PROP(StaggeredModel, GlobalFluxVariablesCache, Dumux::CCGlobalFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
+
+//! The local jacobian operator
+SET_TYPE_PROP(StaggeredModel, LocalJacobian, Dumux::CCLocalJacobian<TypeTag>);
+
+//! Assembler for the global jacobian matrix
+SET_TYPE_PROP(StaggeredModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>);
+
+//! The stencil container
+SET_TYPE_PROP(StaggeredModel, StencilsVector, Dumux::CCStencilsVector<TypeTag>);
+
+//! The local flux variables cache vector class
+SET_TYPE_PROP(StaggeredModel, ElementFluxVariablesCache, Dumux::CCElementFluxVariablesCache<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache)>);
+
+//! The global previous volume variables vector class
+SET_TYPE_PROP(StaggeredModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGlobalVolumeVariablesCache)>);
+
+//! Set the BaseLocalResidual to CCLocalResidual
+SET_TYPE_PROP(StaggeredModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>);
+
+//! indicate that this is no box discretization
+SET_BOOL_PROP(StaggeredModel, ImplicitIsBox, false);
+
+} // namespace Properties
+
+} // namespace Dumux
+
+#endif
diff --git a/test/discretization/CMakeLists.txt b/test/discretization/CMakeLists.txt
index a2f10ebf0a..0b5fb54664 100644
--- a/test/discretization/CMakeLists.txt
+++ b/test/discretization/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory("cellcentered")
+add_subdirectory("staggered")
diff --git a/test/discretization/staggered/CMakeLists.txt b/test/discretization/staggered/CMakeLists.txt
new file mode 100644
index 0000000000..806563abb4
--- /dev/null
+++ b/test/discretization/staggered/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_dumux_test(test_staggered test_staggered test_staggered.cc
+ ${CMAKE_CURRENT_BINARY_DIR}/test_staggered)
diff --git a/test/discretization/staggered/test_staggered.cc b/test/discretization/staggered/test_staggered.cc
new file mode 100644
index 0000000000..1c3a0f9217
--- /dev/null
+++ b/test/discretization/staggered/test_staggered.cc
@@ -0,0 +1,149 @@
+// -*- 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 Test for finite volume element geometry, sub control volume, and sub
+ control volume faces
+ */
+#include <config.h>
+
+#include <iostream>
+#include <utility>
+
+#include <dune/common/test/iteratortest.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dumux/implicit/staggered/properties.hh>
+#include <dumux/discretization/staggered/globalfvgeometry.hh>
+#include <dumux/discretization/staggered/fvelementgeometry.hh>
+#include <dumux/discretization/staggered/subcontrolvolume.hh>
+#include <dumux/discretization/staggered/subcontrolvolumeface.hh>
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class MockProblem
+{
+ using ElementMapper = typename GET_PROP_TYPE(TypeTag, DofMapper);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+public:
+ MockProblem(const GridView& gridView) : mapper_(gridView) {}
+
+ const ElementMapper& elementMapper() const
+ { return mapper_; }
+private:
+ ElementMapper mapper_;
+};
+
+namespace Properties
+{
+NEW_TYPE_TAG(TestFVGeometry, INHERITS_FROM(StaggeredModel));
+
+SET_TYPE_PROP(TestFVGeometry, Grid, Dune::YaspGrid<2>);
+
+SET_TYPE_PROP(TestFVGeometry, Problem, Dumux::MockProblem<TypeTag>);
+
+SET_BOOL_PROP(TestFVGeometry, EnableGlobalFVGeometryCache, true);
+}
+
+}
+
+template<class T>
+class NoopFunctor {
+public:
+ NoopFunctor() {}
+ void operator()(const T& t){}
+};
+
+int main (int argc, char *argv[]) try
+{
+ // maybe initialize mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ std::cout << "Checking the FVGeometries, SCVs and SCV faces" << std::endl;
+
+ // aliases
+ using TypeTag = TTAG(TestFVGeometry);
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using GridView = typename Grid::LeafGridView;
+
+ constexpr int dim = GridView::dimension;
+ constexpr int dimworld = GridView::dimensionworld;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GlobalFVGeometry = typename GET_PROP_TYPE(TypeTag, GlobalFVGeometry);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+
+ // make a grid
+ GlobalPosition lower(0.0);
+ GlobalPosition upper(1.0);
+ std::array<unsigned int, dim> els{{2, 2}};
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+ auto leafGridView = grid->leafGridView();
+
+ Problem problem(leafGridView);
+
+ GlobalFVGeometry global(leafGridView);
+ global.update(problem);
+
+ // iterate over elements. For every element get fv geometry and loop over scvs and scvfaces
+ for (const auto& element : elements(leafGridView))
+ {
+ auto eIdx = problem.elementMapper().index(element);
+ std::cout << std::endl << "Checking fvGeometry of element " << eIdx << std::endl;
+ auto fvGeometry = localView(global);
+ fvGeometry.bind(element);
+
+ auto range = scvs(fvGeometry);
+ NoopFunctor<SubControlVolume> op;
+ if(0 != testForwardIterator(range.begin(), range.end(), op))
+ DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ std::cout << "-- scv " << scv.index() << " center at: " << scv.center() << std::endl;
+ }
+
+ auto range2 = scvfs(fvGeometry);
+ NoopFunctor<SubControlVolumeFace> op2;
+ if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
+ DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal();
+ if (scvf.boundary()) std::cout << " (on boundary).";
+ std::cout << std::endl;
+ }
+ }
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (Dune::Exception e) {
+
+ std::cout << e << std::endl;
+ return 1;
+}
--
GitLab