diff --git a/dumux/implicit/fvelementgeometry.hh b/dumux/implicit/fvelementgeometry.hh new file mode 100644 index 0000000000000000000000000000000000000000..777692987a907044372b52bede41b4372f8b19b8 --- /dev/null +++ b/dumux/implicit/fvelementgeometry.hh @@ -0,0 +1,168 @@ +// -*- 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_FV_ELEMENTGEOMETRY_HH +#define DUMUX_FV_ELEMENTGEOMETRY_HH + +#include <dune/common/iteratorrange.hh> +#include <dune/common/iteratorfacades.hh> + +#include <dumux/implicit/subcontrolvolume.hh> +#include <dumux/implicit/subcontrolvolumeface.hh> + +namespace Dumux +{ +namespace Properties +{ +NEW_PROP_TAG(SubControlVolume); +NEW_PROP_TAG(SubControlVolumeFace); +NEW_PROP_TAG(FVElementGeometry); +NEW_PROP_TAG(FVElementGeometryVector); +} + +/*! + * \ingroup ImplcititModel + * \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: + // Constructor with vectors + FVElementGeometry(const FVElementGeometryVector& fvGeometryVector, + const std::vector<IndexType>& scvIndices, + const std::vector<IndexType>& scvfIndices) + : fvGeometryVector_(fvGeometryVector), scvIndices_(scvIndices), scvfIndices_(scvfIndices) + {} + + // iterator range for sub control volumes + inline Dune::IteratorRange<ScvIterator> scvs() + { + return Dune::IteratorRange<ScvIterator>(ScvIterator(scvIndices_.begin(), fvGeometryVector_), + ScvIterator(scvIndices_.end(), fvGeometryVector_)); + } + + // iterator range for sub control volume faces + inline Dune::IteratorRange<ScvfIterator> scvfs() + { + return Dune::IteratorRange<ScvfIterator>(ScvfIterator(scvfIndices_.begin(), fvGeometryVector_), + ScvfIterator(scvfIndices_.end(), fvGeometryVector_)); + } + +private: + const FVElementGeometryVector& fvGeometryVector_; + std::vector<IndexType> scvIndices_; + std::vector<IndexType> scvfIndices_; +}; + +} + +#endif diff --git a/dumux/implicit/subcontrolvolume.hh b/dumux/implicit/subcontrolvolume.hh new file mode 100644 index 0000000000000000000000000000000000000000..226e37bc5fba827c20e75fde89c01c7b6215e435 --- /dev/null +++ b/dumux/implicit/subcontrolvolume.hh @@ -0,0 +1,111 @@ +// -*- 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_SUBCONTROLVOLUME_HH +#define DUMUX_SUBCONTROLVOLUME_HH + +#include <dune/common/fvector.hh> + +namespace Dumux +{ + +/*! + * \ingroup ImplicitModel + * \brief Base class for a sub control volume, i.e a part of the control + * volume we are making the balance for. + */ +template<class G, typename I> +class SubControlVolume +{ +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: + /*\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 + */ + SubControlVolume(Geometry geometry, + IndexType elementIndex, + IndexType indexInElement) + : geometry_(std::move(geometry)), + elementIndex_(std::move(elementIndex)), + indexInElement_(std::move(indexInElement)) {} + + // The move contructor + SubControlVolume(SubControlVolume&& other) + : geometry_(std::move(other.geometry_)), + elementIndex_(std::move(other.elementIndex_)), + indexInElement_(std::move(other.indexInElement_)) {} + + // The copy contructor + SubControlVolume(const SubControlVolume& other) + : geometry_(other.geometry_), + elementIndex_(other.elementIndex_), + indexInElement_(other.indexInElement_) {} + + //! 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_; + } + + IndexType elementIndex() const + { + return elementIndex_; + } + + IndexType indexInElement() const + { + return indexInElement_; + } + +private: + Geometry geometry_; + IndexType elementIndex_; + IndexType indexInElement_; +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/subcontrolvolumeface.hh b/dumux/implicit/subcontrolvolumeface.hh new file mode 100644 index 0000000000000000000000000000000000000000..2a92474ce4a78bfbc6e17c8fa7553e1f8469cf48 --- /dev/null +++ b/dumux/implicit/subcontrolvolumeface.hh @@ -0,0 +1,120 @@ +// -*- 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_SUBCONTROLVOLUMEFACE_HH +#define DUMUX_SUBCONTROLVOLUMEFACE_HH + +#include <utility> +#include <dune/common/fvector.hh> + +namespace Dumux +{ + +/*! + * \ingroup ImplicitModel + * \brief Base class for a sub control volume face, i.e a part of the boundary + * of a sub control volume we computing a flux on. + */ +template<class Geometry, typename IndexType> +class SubControlVolumeFace +{ + using Scalar = typename Geometry::ctype; + enum { dimworld = Geometry::coorddimension }; + using GlobalPosition = Dune::FieldVector<Scalar, dimworld>; + +public: + SubControlVolumeFace(const Geometry& geometry, + const GlobalPosition& unitOuterNormal, + const std::vector<IndexType>& scvIndices, + const std::vector<IndexType>& volVarsIndices, + bool boundary = false) + : geometry_(geometry), + unitOuterNormal_(unitOuterNormal), + scvIndices_(scvIndices), + volVarsIndices_(volVarsIndices), + boundary_(boundary) {} + + //! The center of the sub control volume face + GlobalPosition center() const + { + return geometry_.center(); + } + + //! The area of the sub control volume face + Scalar area() const + { + return geometry_.volume(); + } + + //! The geometry of the sub control volume face + const Geometry& geometry() const + { + return geometry_; + } + + //! 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 inside sub control volume for upwinding + IndexType insideVolVarsIdx() const + { + return volVarsIndices_[0]; + } + + //! index of the outside sub control volume for spatial param evaluation + // This results in undefined behaviour if boundary is false + IndexType outsideScvIdx() const + { + return scvIndices_[1]; + } + + //! index of the outside sub control volume for upwinding + // This results in undefined behaviour if boundary is false + IndexType outsideVolVarsIdx() const + { + return volVarsIndices_[1]; + } + +private: + Geometry geometry_; + GlobalPosition unitOuterNormal_; + std::vector<IndexType> scvIndices_, volVarsIndices_; + bool boundary_; +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/tpfa/fvelementgeometryvector.hh b/dumux/implicit/tpfa/fvelementgeometryvector.hh new file mode 100644 index 0000000000000000000000000000000000000000..c8edd5c3133be020653b2f3d821101a02839f464 --- /dev/null +++ b/dumux/implicit/tpfa/fvelementgeometryvector.hh @@ -0,0 +1,171 @@ +// -*- 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 cell-centered TPFA models + * This builds up the sub control volumes and sub control volume faces + * for each element. + */ +#ifndef DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH +#define DUMUX_IMPLICIT_TPFA_FV_GEOMETRY_VECTOR_HH + +#include <dumux/implicit/subcontrolvolume.hh> +#include <dumux/implicit/subcontrolvolumeface.hh> +#include <dumux/implicit/fvelementgeometry.hh> + +namespace Dumux +{ +namespace Properties +{ +NEW_PROP_TAG(SubControlVolume); +NEW_PROP_TAG(SubControlVolumeFace); +NEW_PROP_TAG(FVElementGeometry); +NEW_PROP_TAG(FVElementGeometryVector); +} + +//! An index to element map +template <class GridView> +class ElementMap + : public std::vector<typename GridView::Traits::Grid::template Codim<0>::EntitySeed> +{ + using Grid = typename GridView::Traits::Grid; + using Element = typename GridView::template Codim<0>::Entity; + using IndexType = typename GridView::IndexSet::IndexType; +public: + ElementMap(const GridView& gridView_) : grid_(gridView_.grid()) {} + + Element element(IndexType eIdx) + { return grid_.entity((*this)[eIdx]); } + +private: + const Grid& grid_; +}; + +/*! + * \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> +class TpfaFVElementGeometryVector +{ + 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; + +public: + //! Constructor + TpfaFVElementGeometryVector(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. + */ + FVElementGeometry fvGeometry(const Element& element) const + { + return fvGeometry(gridView_.indexSet().index(element)); + } + + FVElementGeometry fvGeometry(IndexType eIdx) const + { + return fvGeometries_[eIdx]; + } + + //! 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]; + } + + // Get an element from a sub control volume contained in it + Element element(const SubControlVolume& scv) const + { return elementMap_.element(scv.elementIndex()); } + + //! update all fvElementGeometries (do this again after grid adaption) + void update() + { + scvs_.clear(); + scvfs_.clear(); + fvGeometries_.clear(); + elementMap_.clear(); + + // Build the SCV and SCV faces + IndexType scvfIdx = 0; + elementMap_.resize(gridView_.size(0)); + scvs_.resize(gridView_.size(0)); + for (const auto& element : elements(gridView_)) + { + auto eIdx = gridView_.indexSet().index(element); + scvs_[eIdx] = std::make_shared<SubControlVolume>(element.geometry(), eIdx, 0); + + // fill the element map with seeds + elementMap_[eIdx] = element.seed(); + + // the element-wise index sets for finite volume geometry + std::vector<IndexType> scvfsIndexSet; + for (const auto& intersection : intersections(gridView_, element)) + { + if (!intersection.boundary()) + { + auto nIdx = gridView_.indexSet().index(intersection.outside()); + scvfs_.push_back(std::make_shared<SubControlVolumeFace>(intersection.geometry(), + intersection.centerUnitOuterNormal(), + std::vector<IndexType>({eIdx, nIdx}), + std::vector<IndexType>({eIdx, nIdx}), + false)); + scvfsIndexSet.push_back(scvfIdx++); + } + else + { + scvfs_.push_back(std::make_shared<SubControlVolumeFace>(intersection.geometry(), + intersection.centerUnitOuterNormal(), + std::vector<IndexType>({eIdx}), + std::vector<IndexType>({eIdx}), + true)); + scvfsIndexSet.push_back(scvfIdx++); + } + } + + // Compute the finite volume element geometries + fvGeometries_.push_back(FVElementGeometry(*this, {eIdx}, scvfsIndexSet)); + } + } + +private: + GridView gridView_; + Dumux::ElementMap<GridView> elementMap_; + std::vector<std::shared_ptr<SubControlVolume>> scvs_; + std::vector<std::shared_ptr<SubControlVolumeFace>> scvfs_; + std::vector<FVElementGeometry> fvGeometries_; +}; + +} // end namespace + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ff547283b232ac8bc0087827cd1b24d50cf1c208..f22a624f481378fa39c4f29a617197418ced7db3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -5,3 +5,4 @@ add_subdirectory("io") add_subdirectory("material") add_subdirectory("multidomain") add_subdirectory("porousmediumflow") +add_subdirectory("implicit") diff --git a/test/implicit/CMakeLists.txt b/test/implicit/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..33078f426d90f1b85dd44ac8a818e337388538d3 --- /dev/null +++ b/test/implicit/CMakeLists.txt @@ -0,0 +1,3 @@ +add_dumux_test(test_fvelementgeometry test_fvelementgeometry test_fvelementgeometry.cc + ${CMAKE_CURRENT_BINARY_DIR}/test_fvelementgeometry) +set(CMAKE_BUILD_TYPE Debug) diff --git a/test/implicit/test_fvelementgeometry.cc b/test/implicit/test_fvelementgeometry.cc new file mode 100644 index 0000000000000000000000000000000000000000..54c1ca1964d3af8d9be339d3efa901f91f64616d --- /dev/null +++ b/test/implicit/test_fvelementgeometry.cc @@ -0,0 +1,154 @@ +// -*- 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 <dumux/common/basicproperties.hh> +#include <dumux/implicit/tpfa/fvelementgeometryvector.hh> +#include <dumux/implicit/fvelementgeometry.hh> +#include <dumux/implicit/subcontrolvolume.hh> +#include <dumux/implicit/subcontrolvolumeface.hh> + +namespace Dumux +{ +namespace Properties +{ +NEW_PROP_TAG(SubControlVolume); +NEW_PROP_TAG(SubControlVolumeFace); +NEW_PROP_TAG(FVElementGeometry); +NEW_PROP_TAG(FVElementGeometryVector); + +NEW_TYPE_TAG(TestFVGeometry, INHERITS_FROM(NumericModel)); + +SET_TYPE_PROP(TestFVGeometry, Grid, Dune::YaspGrid<2>); + +SET_TYPE_PROP(TestFVGeometry, GridView, typename GET_PROP_TYPE(TypeTag, Grid)::LeafGridView); + +SET_PROP(TestFVGeometry, 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::SubControlVolume<ScvGeometry, IndexType> type; +}; + +SET_PROP(TestFVGeometry, 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::SubControlVolumeFace<ScvfGeometry, IndexType> type; + +}; + +SET_TYPE_PROP(TestFVGeometry, FVElementGeometry, FVElementGeometry<TypeTag>); + +SET_TYPE_PROP(TestFVGeometry, FVElementGeometryVector, TpfaFVElementGeometryVector<TypeTag>); + +} +} + +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 FVElementGeometryVector = typename GET_PROP_TYPE(TypeTag, FVElementGeometryVector); + + // 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(); + + FVElementGeometryVector fvGeometries(leafGridView); + fvGeometries.update(); + + // iterate over elements. For every element get fv geometry and loop over scvs and scvfaces + for (const auto& element : elements(leafGridView)) + { + std::cout << std::endl << "Checking fvGeometry of element " << leafGridView.indexSet().index(element) << std::endl; + auto fvGeometry = fvGeometries.fvGeometry(element); + + auto range = fvGeometry.scvs(); + 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 : fvGeometry.scvs()) + { + std::cout << "-- scv center at: " << scv.center() << std::endl; + } + + auto range2 = fvGeometry.scvfs(); + 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 : fvGeometry.scvfs()) + { + std::cout << "-- scvf center at: " << scvf.center() << std::endl; + } + } +} +// ////////////////////////////////// +// Error handler +// ///////////////////////////////// +catch (Dune::Exception e) { + + std::cout << e << std::endl; + return 1; +}