diff --git a/CHANGELOG.md b/CHANGELOG.md
index d0c6549ae0c3aec0c0f7794d42d284c9bf96cf25..59e37b6510e065682a50daddeccef1d06695e442 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,8 @@ Differences Between DuMuX 3.2 and DuMuX 3.1
### Improvements and Enhancements
+- __Radially symmetric problems__: We now have support for radially symmetric problems (disc, ball, toroid). The support comes in form of wrappers for sub control volumes and faces that overload the respective `volume()` and `area()` function turning a 1d or 2d problem into a 2d or 3d radially symmetric problem.
+
### Immediate interface changes not allowing/requiring a deprecation period
### Deprecated properties, to be removed after 3.2:
diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt
index c12d297aa602f617267f96848505a36960b75033..b2a12576239b166f62735478ce043d1f0d021718 100644
--- a/dumux/discretization/CMakeLists.txt
+++ b/dumux/discretization/CMakeLists.txt
@@ -20,6 +20,9 @@ fvgridvariables.hh
fvproperties.hh
localview.hh
method.hh
+rotationpolicy.hh
+rotationsymmetricscv.hh
+rotationsymmetricscvf.hh
scvandscvfiterators.hh
staggered.hh
subcontrolvolumebase.hh
diff --git a/dumux/discretization/rotationpolicy.hh b/dumux/discretization/rotationpolicy.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9e762ef6edb22fc3eb15f2242637ffd63cb52a98
--- /dev/null
+++ b/dumux/discretization/rotationpolicy.hh
@@ -0,0 +1,41 @@
+// -*- 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 3 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
+ * \ingroup Discretization
+ * \brief Rotation policy for defining rotational symmetric grid geometries
+ */
+#ifndef DUMUX_DISCRETIZATION_ROTATION_POLICY_HH
+#define DUMUX_DISCRETIZATION_ROTATION_POLICY_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \brief Rotation policies
+ * - disc (or annulus): rotate a segment around an axis perpendicular to the line through the segment (we choose )
+ * - ball (or shell): rotate a segment around a point on the line through the segment
+ * - toroid: rotate a polygon around one of its axis (we rotate about the second axis)
+ */
+enum class RotationPolicy
+{ disc, ball, toroid };
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/rotationsymmetricscv.hh b/dumux/discretization/rotationsymmetricscv.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0870dc6445d0691687729e7cd7864b87dd6eb7c0
--- /dev/null
+++ b/dumux/discretization/rotationsymmetricscv.hh
@@ -0,0 +1,123 @@
+// -*- 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 3 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
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume rotation symmetric
+ */
+#ifndef DUMUX_DISCRETIZATION_ROTATION_SYMMETRIC_SUBCONTROLVOLUME_HH
+#define DUMUX_DISCRETIZATION_ROTATION_SYMMETRIC_SUBCONTROLVOLUME_HH
+
+#include <cmath>
+#include <dumux/discretization/rotationpolicy.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume rotation symmetric
+ * \tparam SubControlVolume The wrapped scv type
+ * \tparam rotationPolicy the rotation policy (see enum RotationPolicy)
+ */
+template<class SubControlVolume, RotationPolicy rotationPolicy>
+class RotationSymmetricSubControlVolume;
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume rotation symmetric
+ * \tparam SubControlVolume The wrapped scv type
+ * \note Specialization for the 'disc' policy (1d grid --> 2d disc)
+ */
+template<class SubControlVolume>
+class RotationSymmetricSubControlVolume<SubControlVolume, RotationPolicy::disc>
+: public SubControlVolume
+{
+ using Scalar = typename SubControlVolume::Traits::Scalar;
+ static_assert(SubControlVolume::Traits::Geometry::mydimension == 1, "Rotation symmetric scv with disc policy only works with 1d grids!");
+ static_assert(SubControlVolume::Traits::Geometry::coorddimension == 1, "Rotation symmetric scv with disc policy only works with 1d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolume::SubControlVolume;
+
+ //! The volume of the sub control volume (difference between two disks)
+ Scalar volume() const
+ {
+ using std::abs;
+ const auto radius0 = this->corner(0)[0];
+ const auto radius1 = this->corner(1)[0];
+ return M_PI*abs(radius1*radius1 - radius0*radius0);
+ }
+};
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume rotation symmetric
+ * \tparam SubControlVolume The wrapped scv type
+ * \note Specialization for the 'ball' policy (1d grid --> 3d ball)
+ */
+template<class SubControlVolume>
+class RotationSymmetricSubControlVolume<SubControlVolume, RotationPolicy::ball>
+: public SubControlVolume
+{
+ using Scalar = typename SubControlVolume::Traits::Scalar;
+ static_assert(SubControlVolume::Traits::Geometry::mydimension == 1, "Rotation symmetric scv with ball policy only works with 1d grids!");
+ static_assert(SubControlVolume::Traits::Geometry::coorddimension == 1, "Rotation symmetric scv with ball policy only works with 1d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolume::SubControlVolume;
+
+ //! The volume of the sub control volume (difference between two balls)
+ Scalar volume() const
+ {
+ using std::abs;
+ const auto radius0 = this->corner(0)[0];
+ const auto radius1 = this->corner(1)[0];
+ return 4.0/3.0*M_PI*abs(radius1*radius1*radius1 - radius0*radius0*radius0);
+ }
+};
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume rotation symmetric
+ * \tparam SubControlVolume The wrapped scv type
+ * \note Specialization for the 'toroid' policy (2d grid --> 3d toroid)
+ * \note We rotate about the second axis
+ */
+template<class SubControlVolume>
+class RotationSymmetricSubControlVolume<SubControlVolume, RotationPolicy::toroid>
+: public SubControlVolume
+{
+ using Scalar = typename SubControlVolume::Traits::Scalar;
+ static_assert(SubControlVolume::Traits::Geometry::mydimension == 2, "Rotation symmetric scv with toroid policy only works with 2d grids!");
+ static_assert(SubControlVolume::Traits::Geometry::coorddimension == 2, "Rotation symmetric scv with toroid policy only works with 2d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolume::SubControlVolume;
+
+ //! The volume of the sub control volume (Guldinus theorem)
+ Scalar volume() const
+ {
+ const auto radius = this->center()[1];
+ return SubControlVolume::volume()*2.0*M_PI*radius;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/rotationsymmetricscvf.hh b/dumux/discretization/rotationsymmetricscvf.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5075fedf5dc41cdd1baf18fcabf276860c1b3d51
--- /dev/null
+++ b/dumux/discretization/rotationsymmetricscvf.hh
@@ -0,0 +1,119 @@
+// -*- 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 3 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
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume face rotation symmetric
+ */
+#ifndef DUMUX_DISCRETIZATION_ROTATION_SYMMETRIC_SUBCONTROLVOLUMEFACE_HH
+#define DUMUX_DISCRETIZATION_ROTATION_SYMMETRIC_SUBCONTROLVOLUMEFACE_HH
+
+#include <cmath>
+#include <dumux/discretization/rotationpolicy.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume face rotation symmetric
+ * \tparam SubControlVolumeFace The wrapped scvf type
+ * \tparam rotationPolicy the rotation policy (see enum RotationPolicy)
+ */
+template<class SubControlVolumeFace, RotationPolicy rotationPolicy>
+class RotationSymmetricSubControlVolumeFace;
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume face rotation symmetric
+ * \tparam SubControlVolumeFace The wrapped scvf type
+ * \note Specialization for the 'disc' policy (1d grid --> 2d disc)
+ */
+template<class SubControlVolumeFace>
+class RotationSymmetricSubControlVolumeFace<SubControlVolumeFace, RotationPolicy::disc>
+: public SubControlVolumeFace
+{
+ using Scalar = typename SubControlVolumeFace::Traits::Scalar;
+ static_assert(SubControlVolumeFace::Traits::Geometry::mydimension == 0, "Rotation symmetric scvf with disc policy only works with 1d grids!");
+ static_assert(SubControlVolumeFace::Traits::Geometry::coorddimension == 1, "Rotation symmetric scvf with disc policy only works with 1d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolumeFace::SubControlVolumeFace;
+
+ //! The area of the sub control volume face (circumference of circle)
+ Scalar area() const
+ {
+ const auto radius = this->corner(0)[0];
+ return 2.0*M_PI*radius;
+ }
+};
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume face rotation symmetric
+ * \tparam SubControlVolumeFace The wrapped scvf type
+ * \note Specialization for the 'ball' policy (1d grid --> 3d ball)
+ */
+template<class SubControlVolumeFace>
+class RotationSymmetricSubControlVolumeFace<SubControlVolumeFace, RotationPolicy::ball>
+: public SubControlVolumeFace
+{
+ using Scalar = typename SubControlVolumeFace::Traits::Scalar;
+ static_assert(SubControlVolumeFace::Traits::Geometry::mydimension == 0, "Rotation symmetric scvf with ball policy only works with 1d grids!");
+ static_assert(SubControlVolumeFace::Traits::Geometry::coorddimension == 1, "Rotation symmetric scvf with ball policy only works with 1d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolumeFace::SubControlVolumeFace;
+
+ //! The area of the sub control volume face (surface of sphere)
+ Scalar area() const
+ {
+ const auto radius = this->corner(0)[0];
+ return 4.0*M_PI*radius*radius;
+ }
+};
+
+/*!
+ * \ingroup Discretization
+ * \brief Wrapper to make a sub control volume face rotation symmetric
+ * \tparam SubControlVolumeFace The wrapped scvf type
+ * \note Specialization for the 'toroid' policy (2d grid --> 3d toroid)
+ * \note We rotate about the second axis
+ */
+template<class SubControlVolumeFace>
+class RotationSymmetricSubControlVolumeFace<SubControlVolumeFace, RotationPolicy::toroid>
+: public SubControlVolumeFace
+{
+ using Scalar = typename SubControlVolumeFace::Traits::Scalar;
+ static_assert(SubControlVolumeFace::Traits::Geometry::mydimension == 1, "Rotation symmetric scvf with toroid policy only works with 2d grids!");
+ static_assert(SubControlVolumeFace::Traits::Geometry::coorddimension == 2, "Rotation symmetric scvf with toroid policy only works with 2d grids!");
+public:
+ //! Parent type constructor
+ using SubControlVolumeFace::SubControlVolumeFace;
+
+ //! The area of the sub control volume face (Guldinus theorem)
+ Scalar area() const
+ {
+ const auto radius = this->center()[1];
+ return SubControlVolumeFace::area()*2.0*M_PI*radius;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/discretization/CMakeLists.txt b/test/discretization/CMakeLists.txt
index f127b603f79c94575badc24a7eb39c9328c88892..25c1bdd297e74a772cd7d985038aa7cde682c038 100644
--- a/test/discretization/CMakeLists.txt
+++ b/test/discretization/CMakeLists.txt
@@ -2,3 +2,4 @@ add_subdirectory(cellcentered)
add_subdirectory(staggered)
add_subdirectory(box)
add_subdirectory(projection)
+add_subdirectory(rotationsymmetry)
diff --git a/test/discretization/rotationsymmetry/CMakeLists.txt b/test/discretization/rotationsymmetry/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a8237b62853ccb9adcd4b836bc631644bef27442
--- /dev/null
+++ b/test/discretization/rotationsymmetry/CMakeLists.txt
@@ -0,0 +1,2 @@
+dune_add_test(SOURCES test_rotationsymmetric_gridgeometry.cc
+ LABELS unit discretization)
diff --git a/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc
new file mode 100644
index 0000000000000000000000000000000000000000..15072fb5227dfca3f54eb95dc4881cb5d35ca793
--- /dev/null
+++ b/test/discretization/rotationsymmetry/test_rotationsymmetric_gridgeometry.cc
@@ -0,0 +1,208 @@
+// -*- 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 3 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 rotation symmetric grid geometry
+ */
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/float_cmp.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
+#include <dumux/discretization/rotationpolicy.hh>
+#include <dumux/discretization/rotationsymmetricscv.hh>
+#include <dumux/discretization/rotationsymmetricscvf.hh>
+
+namespace Dumux {
+
+template<class BaseTraits, RotationPolicy rotPolicy>
+struct GGRotSymTraits : public BaseTraits
+{
+ using SubControlVolume = RotationSymmetricSubControlVolume<typename BaseTraits::SubControlVolume, rotPolicy>;
+ using SubControlVolumeFace = RotationSymmetricSubControlVolumeFace<typename BaseTraits::SubControlVolumeFace, rotPolicy>;
+};
+
+template<class GG>
+void runTest(const GG& gg, const double refVolume, const double refSurface)
+{
+ double volume = 0.0;
+ double surface = 0.0;
+ for (const auto& element : elements(gg.gridView()))
+ {
+ auto fvGeometry = localView(gg);
+ fvGeometry.bind(element);
+
+ for (const auto& scv : scvs(fvGeometry))
+ volume += scv.volume();
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (scvf.boundary())
+ surface += scvf.area();
+ }
+
+ // compare to reference
+ if (!Dune::FloatCmp::eq(refVolume, volume))
+ DUNE_THROW(Dune::Exception, "Volume not correct! Reference: " << refVolume << ", computed: " << volume);
+ if (!Dune::FloatCmp::eq(refSurface, surface))
+ DUNE_THROW(Dune::Exception, "Surface not correct! Reference: " << refSurface << ", computed: " << surface);
+}
+
+} // end namespace Dumux
+
+int main (int argc, char *argv[]) try
+{
+ using namespace Dumux;
+
+ // maybe initialize mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ // test the disc policy
+ {
+ using Grid = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>;
+ using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::disc>;
+ using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>;
+ using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition;
+
+ // make a grid
+ const double innerRadius = 0.1;
+ const double outerRadius = 1.0;
+ GlobalPosition lower(innerRadius);
+ GlobalPosition upper(outerRadius);
+ std::array<unsigned int, Grid::dimension> els{{10}};
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+
+ // obtain leaf and make GridGeometry
+ auto leafGridView = grid->leafGridView();
+ GridGeometry gg(leafGridView);
+ gg.update();
+
+ // compute the annulus area and the surface
+ const double refVolume = M_PI*(outerRadius*outerRadius - innerRadius*innerRadius);
+ const double refSurface = 2.0*M_PI*(innerRadius + outerRadius);
+ runTest(gg, refVolume, refSurface);
+
+ std::cout << "Successfully tested disc (annulus) policy." << std::endl;
+
+ } // end disc policy
+
+ // test the ball policy
+ {
+ using Grid = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>;
+ using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::ball>;
+ using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>;
+ using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition;
+
+ // make a grid
+ const double innerRadius = 0.1;
+ const double outerRadius = 1.0;
+ GlobalPosition lower(innerRadius);
+ GlobalPosition upper(outerRadius);
+ std::array<unsigned int, Grid::dimension> els{{10}};
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+
+ // obtain leaf and make GridGeometry
+ auto leafGridView = grid->leafGridView();
+ GridGeometry gg(leafGridView);
+ gg.update();
+
+ // compute the annulus area and the surface
+ const double refVolume = 4.0/3.0*M_PI*(outerRadius*outerRadius*outerRadius - innerRadius*innerRadius*innerRadius);
+ const double refSurface = 4.0*M_PI*(innerRadius*innerRadius + outerRadius*outerRadius);
+ runTest(gg, refVolume, refSurface);
+
+ std::cout << "Successfully tested ball (shell) policy." << std::endl;
+
+ } // end ball policy
+
+ // test the toroid policy
+ {
+ using Grid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>;
+ using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::toroid>;
+ using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>;
+ using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition;
+
+ // make a grid
+ const double innerRadius = 0.1;
+ const double outerRadius = 1.0;
+ GlobalPosition lower({innerRadius, innerRadius});
+ GlobalPosition upper({outerRadius, outerRadius});
+ std::array<unsigned int, Grid::dimension> els({20, 20});
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+
+ // obtain leaf and make GridGeometry
+ auto leafGridView = grid->leafGridView();
+ GridGeometry gg(leafGridView);
+ gg.update();
+
+ // compute the annulus area and the surface
+ const auto centroidRadius = 0.5*(innerRadius + outerRadius);
+ const auto side = (outerRadius - innerRadius);
+ const double refVolume = side*side*2.0*M_PI*centroidRadius;
+ const double refSurface = 4.0*side*2.0*M_PI*centroidRadius;
+ runTest(gg, refVolume, refSurface);
+
+ std::cout << "Successfully tested toroid policy." << std::endl;
+
+ } // end toroid policy
+
+ // test the toroid policy for perfect cylinder
+ {
+ using Grid = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>;
+ using GGTraits = GGRotSymTraits<CCTpfaDefaultGridGeometryTraits<typename Grid::LeafGridView>, RotationPolicy::toroid>;
+ using GridGeometry = CCTpfaFVGridGeometry<typename Grid::LeafGridView, /*caching=*/false, GGTraits>;
+ using GlobalPosition = typename GridGeometry::SubControlVolume::GlobalPosition;
+
+ // make a grid
+ const double innerRadius = 0.0;
+ const double outerRadius = 1.0;
+ GlobalPosition lower({innerRadius, innerRadius});
+ GlobalPosition upper({outerRadius, outerRadius});
+ std::array<unsigned int, Grid::dimension> els({20, 20});
+ std::shared_ptr<Grid> grid = Dune::StructuredGridFactory<Grid>::createCubeGrid(lower, upper, els);
+
+ // obtain leaf and make GridGeometry
+ auto leafGridView = grid->leafGridView();
+ GridGeometry gg(leafGridView);
+ gg.update();
+
+ // compute the annulus area and the surface
+ const double refVolume = outerRadius*M_PI*outerRadius*outerRadius;
+ const double refSurface = 2.0*M_PI*outerRadius + 2.0*M_PI*outerRadius*outerRadius;
+ runTest(gg, refVolume, refSurface);
+
+ std::cout << "Successfully tested toroid policy for perfect cylinder." << std::endl;
+
+ } // end toroid policy
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (Dune::Exception &e) {
+
+ std::cout << e << std::endl;
+ return 1;
+}