From bc9d789924ab789e08831716348b29de74e36787 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:46:38 +0200
Subject: [PATCH 01/10] [cleanup] Remove unnecessary includes

---
 dumux/freeflow/navierstokes/staggered/fluxoversurface.hh        | 1 -
 .../2p/sequential/diffusion/cellcentered/velocityadaptive.hh    | 2 --
 dumux/porousmediumflow/velocityoutput.hh                        | 2 --
 test/io/gridmanager/gmshboundaryflagtest.hh                     | 1 -
 test/multidomain/embedded/1d3d/1p2c_richards2c/problem_soil.hh  | 1 -
 test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh          | 1 -
 test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh      | 1 -
 test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh          | 1 -
 test/multidomain/facet/test_gridmanager.cc                      | 1 -
 9 files changed, 11 deletions(-)

diff --git a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
index 30d3e638ec..79d7474006 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
@@ -33,7 +33,6 @@
 #include <dune/common/fvector.hh>
 #include <dune/geometry/type.hh>
 #include <dune/geometry/multilineargeometry.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/geometry/makegeometry.hh>
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocityadaptive.hh
index 0e06a573cf..e984fcd6e2 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocityadaptive.hh
@@ -59,8 +59,6 @@ class FVVelocity2PAdaptive: public FVVelocity2P<TypeTag>
         dim = GridView::dimension, dimWorld = GridView::dimensionworld
     };
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     enum
     {
         pw = Indices::pressureW,
diff --git a/dumux/porousmediumflow/velocityoutput.hh b/dumux/porousmediumflow/velocityoutput.hh
index aec441d395..8f2a77ebc5 100644
--- a/dumux/porousmediumflow/velocityoutput.hh
+++ b/dumux/porousmediumflow/velocityoutput.hh
@@ -27,7 +27,6 @@
 
 #include <memory>
 #include <dune/common/float_cmp.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/parameters.hh>
 #include <dumux/io/velocityoutput.hh>
@@ -64,7 +63,6 @@ class PorousMediumFlowVelocityOutput : public VelocityOutput<GridVariables>
     static constexpr int dofCodim = isBox ? dim : 0;
 
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    using ReferenceElements = Dune::ReferenceElements<typename GridView::ctype, dim>;
 
     using Problem = typename GridVolumeVariables::Problem;
     using BoundaryTypes = typename Problem::Traits::BoundaryTypes;
diff --git a/test/io/gridmanager/gmshboundaryflagtest.hh b/test/io/gridmanager/gmshboundaryflagtest.hh
index 4d242a469b..4f67a03803 100644
--- a/test/io/gridmanager/gmshboundaryflagtest.hh
+++ b/test/io/gridmanager/gmshboundaryflagtest.hh
@@ -40,7 +40,6 @@ class GmshBoundaryFlagTest
     using Scalar = double;
     static const int dim = Grid::dimension;
     using GridManager = typename Dumux::GridManager<Grid>;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
 
 public:
 
diff --git a/test/multidomain/embedded/1d3d/1p2c_richards2c/problem_soil.hh b/test/multidomain/embedded/1d3d/1p2c_richards2c/problem_soil.hh
index ccde53d420..322a7dda40 100644
--- a/test/multidomain/embedded/1d3d/1p2c_richards2c/problem_soil.hh
+++ b/test/multidomain/embedded/1d3d/1p2c_richards2c/problem_soil.hh
@@ -27,7 +27,6 @@
 #define DUMUX_TISSUE_PROBLEM_HH
 
 #include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/grid/yaspgrid.hh>
 #include <dune/grid/uggrid.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
diff --git a/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh b/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh
index 58d76aeb0b..4202ca8346 100644
--- a/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh
+++ b/test/multidomain/embedded/1d3d/1p_1p/problem_tissue.hh
@@ -29,7 +29,6 @@
 #include <dune/grid/yaspgrid.hh>
 
 #include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/common/math.hh>
diff --git a/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh b/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
index 9d0fc7f9bd..4257a04421 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
+++ b/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
@@ -28,7 +28,6 @@
 
 #include <dune/grid/yaspgrid.hh>
 #include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/common/math.hh>
diff --git a/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh b/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh
index 1365c670ee..4e18bf1a61 100644
--- a/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh
+++ b/test/multidomain/embedded/2d3d/1p_1p/problem_matrix.hh
@@ -26,7 +26,6 @@
 #define DUMUX_MATRIX_ROBLEM_HH
 
 #include <dune/geometry/quadraturerules.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/grid/yaspgrid.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
diff --git a/test/multidomain/facet/test_gridmanager.cc b/test/multidomain/facet/test_gridmanager.cc
index 8a14d9725c..ce40d13cae 100644
--- a/test/multidomain/facet/test_gridmanager.cc
+++ b/test/multidomain/facet/test_gridmanager.cc
@@ -30,7 +30,6 @@
 
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/common/mcmgmapper.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/alugrid/grid.hh>
 #include <dune/foamgrid/foamgrid.hh>
 
-- 
GitLab


From 65ca72909b353c84a749f2a0021f145c4a7e3bfc Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:58:51 +0200
Subject: [PATCH 02/10] [io][vtkmultiwriter] replace all
 ReferenceElements::general calls with referenceElement

---
 dumux/io/vtkmultiwriter.hh | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/dumux/io/vtkmultiwriter.hh b/dumux/io/vtkmultiwriter.hh
index 573326cdb6..1438db0723 100644
--- a/dumux/io/vtkmultiwriter.hh
+++ b/dumux/io/vtkmultiwriter.hh
@@ -59,7 +59,6 @@ class VtkNestedFunction : public Dune::VTKFunction<GridView>
     enum { dim = GridView::dimension };
     using ctype = typename GridView::ctype;
     using Element = typename GridView::template Codim<0>::Entity;
-    using ReferenceElements = Dune::ReferenceElements<ctype, dim>;
 
 public:
     VtkNestedFunction(std::string name,
@@ -98,13 +97,11 @@ public:
             // coordinates. This code is based on Dune::P1VTKFunction
             double min=1e100;
             int imin=-1;
-            Dune::GeometryType geomType = element.type();
             int n = element.subEntities(dim);
 
             for (int i=0; i < n; ++i)
             {
-                Dune::FieldVector<ctype,dim> local =
-                    ReferenceElements::general(geomType).position(i,dim);
+                Dune::FieldVector<ctype,dim> local = referenceElement(element).position(i,dim);
                 local -= xi;
                 if (local.infinity_norm()<min)
                 {
-- 
GitLab


From a614e1256a01626c7130cd27b31abd06827c5779 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:47:14 +0200
Subject: [PATCH 03/10] [pmf][velocity] replace all ReferenceElements::general
 calls with referenceElement

---
 dumux/porousmediumflow/velocity.hh | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh
index c8393439e8..f47142fb4d 100644
--- a/dumux/porousmediumflow/velocity.hh
+++ b/dumux/porousmediumflow/velocity.hh
@@ -26,15 +26,15 @@
 #define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
 
 #include <vector>
-#include <dumux/flux/traits.hh>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/float_cmp.hh>
-#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
 
 #include <dumux/common/parameters.hh>
 #include <dumux/discretization/method.hh>
 #include <dumux/discretization/elementsolution.hh>
+#include <dumux/flux/traits.hh>
 
 namespace Dumux {
 
@@ -67,7 +67,6 @@ class PorousMediumFlowVelocity
     static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
 
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    using ReferenceElements = Dune::ReferenceElements<typename GridView::ctype, dim>;
 
     using Problem = typename GridVolumeVariables::Problem;
     using BoundaryTypes = typename Problem::Traits::BoundaryTypes;
@@ -145,8 +144,8 @@ public:
         }
 
         // get the transposed Jacobian of the element mapping
-        const auto referenceElement = ReferenceElements::general(geomType);
-        const auto& localPos = referenceElement.position(0, 0);
+        const auto refElement = referenceElement(geometry);
+        const auto& localPos = refElement.position(0, 0);
         const auto jacobianT2 = geometry.jacobianTransposed(localPos);
 
         if(isBox)
-- 
GitLab


From 43dbcd243177dc9985f2c21af1c46c7514b2ab45 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:43:01 +0200
Subject: [PATCH 04/10] [discretization][box] replace all
 ReferenceElements::general calls with referenceElement

---
 dumux/discretization/box/boxgeometryhelper.hh | 32 +++++++--------
 dumux/discretization/box/fvelementgeometry.hh | 11 ++----
 dumux/discretization/box/fvgridgeometry.hh    | 39 ++++++++-----------
 3 files changed, 35 insertions(+), 47 deletions(-)

diff --git a/dumux/discretization/box/boxgeometryhelper.hh b/dumux/discretization/box/boxgeometryhelper.hh
index c99c5ecbad..be6144326a 100644
--- a/dumux/discretization/box/boxgeometryhelper.hh
+++ b/dumux/discretization/box/boxgeometryhelper.hh
@@ -26,7 +26,6 @@
 #define DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
 
 #include <array>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/math.hh>
 
@@ -128,7 +127,6 @@ class BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
 
     static constexpr auto dim = GridView::dimension;
     static constexpr auto dimWorld = GridView::dimensionworld;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
 
     //! the maximum number of helper points used to construct the geometries
     //! Using a statically sized point array is much faster than dynamic allocation
@@ -140,7 +138,7 @@ public:
     : elementGeometry_(geometry)
     , corners_(geometry.corners())
     {
-        const auto referenceElement = ReferenceElements::general(geometry.type());
+        const auto refElement = referenceElement(geometry);
 
         // the element center
         p_[0] = geometry.center();
@@ -150,8 +148,8 @@ public:
             p_[i+1] = geometry.corner(i);
 
         // face midpoints
-        for (int i = 0; i < referenceElement.size(1); ++i)
-            p_[i+corners_+1] = geometry.global(referenceElement.position(i, 1));
+        for (int i = 0; i < refElement.size(1); ++i)
+            p_[i+corners_+1] = geometry.global(refElement.position(i, 1));
     }
 
     //! Create a vector with the scv corners
@@ -250,9 +248,9 @@ public:
                                              const typename Intersection::Geometry& isGeom,
                                              unsigned int indexInIntersection) const
     {
-        const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
+        const auto refElement = referenceElement(elementGeometry_);
 
-        const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, indexInIntersection, dim);
+        const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, indexInIntersection, dim);
         if (indexInIntersection == 0)
             return ScvfCornerStorage({p_[vIdxLocal+1], isGeom.center()});
         else if (indexInIntersection == 1)
@@ -348,8 +346,6 @@ class BoxGeometryHelper<GridView, 3, ScvType, ScvfType>
 
     static constexpr auto dim = GridView::dimension;
     static constexpr auto dimWorld = GridView::dimensionworld;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
-    using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
 
     //! the maximum number of helper points used to construct the geometries
     //! Using a statically sized point array is much faster than dynamic allocation
@@ -359,7 +355,7 @@ public:
     : elementGeometry_(geometry)
     , corners_(geometry.corners())
     {
-        const auto referenceElement = ReferenceElements::general(geometry.type());
+        const auto refElement = referenceElement(geometry);
 
         // the element center
         p_[0] = geometry.center();
@@ -369,12 +365,12 @@ public:
             p_[i+1] = geometry.corner(i);
 
         // edge midpoints
-        for (int i = 0; i < referenceElement.size(dim-1); ++i)
-            p_[i+corners_+1] = geometry.global(referenceElement.position(i, dim-1));
+        for (int i = 0; i < refElement.size(dim-1); ++i)
+            p_[i+corners_+1] = geometry.global(refElement.position(i, dim-1));
 
         // face midpoints
-        for (int i = 0; i < referenceElement.size(1); ++i)
-            p_[i+corners_+1+referenceElement.size(dim-1)] = geometry.global(referenceElement.position(i, 1));
+        for (int i = 0; i < refElement.size(1); ++i)
+            p_[i+corners_+1+refElement.size(dim-1)] = geometry.global(refElement.position(i, 1));
     }
 
     //! Create a vector with the scv corners
@@ -552,8 +548,8 @@ public:
                                              const typename Intersection::Geometry& geometry,
                                              unsigned int indexInIntersection) const
     {
-        const auto referenceElement = ReferenceElements::general(elementGeometry_.type());
-        const auto faceRefElem = FaceReferenceElements::general(geometry.type());
+        const auto refElement = referenceElement(elementGeometry_);
+        const auto faceRefElem = referenceElement(geometry);
 
         GlobalPosition pi[9];
         auto corners = geometry.corners();
@@ -565,14 +561,14 @@ public:
         const auto idxInInside = is.indexInInside();
         for (int i = 0; i < corners; ++i)
         {
-            const auto vIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim);
+            const auto vIdxLocal = refElement.subEntity(idxInInside, 1, i, dim);
             pi[i+1] = elementGeometry_.corner(vIdxLocal);
         }
 
         // edge midpoints
         for (int i = 0; i < faceRefElem.size(1); ++i)
         {
-            const auto edgeIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim-1);
+            const auto edgeIdxLocal = refElement.subEntity(idxInInside, 1, i, dim-1);
             pi[i+corners+1] = p_[edgeIdxLocal+corners_+1];
         }
 
diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index 9ca3a412bf..1dc89bca1f 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -27,7 +27,6 @@
 #define DUMUX_DISCRETIZATION_BOX_FV_ELEMENT_GEOMETRY_HH
 
 #include <dune/geometry/type.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/common/indextraits.hh>
@@ -58,7 +57,6 @@ class BoxFVElementGeometry<GG, true>
     using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 public:
     //! export the element type
     using Element = typename GridView::template Codim<0>::Entity;
@@ -174,7 +172,6 @@ class BoxFVElementGeometry<GG, false>
     using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
     using GeometryHelper = BoxGeometryHelper<GridView, dim,
                                              typename GG::SubControlVolume,
@@ -283,7 +280,7 @@ private:
         // get the element geometry
         auto elementGeometry = element.geometry();
         elemGeometryType_ = elementGeometry.type();
-        const auto referenceElement = ReferenceElements::general(elemGeometryType_);
+        const auto refElement = referenceElement(elementGeometry);
 
         // get the sub control volume geometries of this element
         GeometryHelper geometryHelper(elementGeometry);
@@ -310,8 +307,8 @@ private:
         for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
         {
             // find the local scv indices this scvf is connected to
-            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                         static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                         static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
             scvfs_[scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
                                                         element,
@@ -332,7 +329,7 @@ private:
                 for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
                 {
                     // find the scv this scvf is connected to
-                    const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
+                    const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
                     std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
 
                     scvfs_.emplace_back(geometryHelper,
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index fd1ecd536c..a1cea412c4 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -28,7 +28,6 @@
 
 #include <unordered_map>
 
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/discretization/method.hh>
@@ -91,8 +90,6 @@ class BoxFVGridGeometry<Scalar, GV, true, Traits>
     static const int dim = GV::dimension;
     static const int dimWorld = GV::dimensionworld;
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
     using GeometryHelper = BoxGeometryHelper<GV, dim,
                                              typename Traits::SubControlVolume,
                                              typename Traits::SubControlVolumeFace>;
@@ -171,7 +168,7 @@ public:
 
             // get the element geometry
             auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // instantiate the geometry helper
             GeometryHelper geometryHelper(elementGeometry);
@@ -194,8 +191,8 @@ public:
             for (; scvfLocalIdx < element.subEntities(dim-1); ++scvfLocalIdx)
             {
                 // find the global and local scv indices this scvf is belonging to
-                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                             static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                             static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
                 scvfs_[eIdx][scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
                                                                   element,
@@ -220,7 +217,7 @@ public:
                     for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
                     {
                         // find the scvs this scvf is belonging to
-                        const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
+                        const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(intersection.indexInInside(), 1, isScvfLocalIdx, dim));
                         std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
 
                         scvfs_[eIdx].emplace_back(geometryHelper,
@@ -238,10 +235,10 @@ public:
                     // add all vertices on the intersection to the set of
                     // boundary vertices
                     const auto fIdx = intersection.indexInInside();
-                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
@@ -254,11 +251,11 @@ public:
 
                     // find the mapped periodic vertex of all vertices on periodic boundaries
                     const auto fIdx = intersection.indexInInside();
-                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                     const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         const auto vPos = elementGeometry.corner(vIdx);
 
@@ -270,10 +267,10 @@ public:
                             if (isOutside.boundary() && isOutside.neighbor())
                             {
                                 const auto fIdxOutside = isOutside.indexInInside();
-                                const auto numFaceVertsOutside = referenceElement.size(fIdxOutside, 1, dim);
+                                const auto numFaceVertsOutside = refElement.size(fIdxOutside, 1, dim);
                                 for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
                                 {
-                                    const auto vIdxOutside = referenceElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
+                                    const auto vIdxOutside = refElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
                                     const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
                                     const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
                                     if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
@@ -363,8 +360,6 @@ class BoxFVGridGeometry<Scalar, GV, false, Traits>
     using Element = typename GV::template Codim<0>::Entity;
     using CoordScalar = typename GV::ctype;
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
 public:
     //! export discretization method
     static constexpr DiscretizationMethod discMethod = DiscretizationMethod::box;
@@ -427,7 +422,7 @@ public:
             numScvf_ += element.subEntities(dim-1);
 
             const auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // store the sub control volume face indices on the domain boundary
             for (const auto& intersection : intersections(this->gridView(), element))
@@ -441,10 +436,10 @@ public:
                     // add all vertices on the intersection to the set of
                     // boundary vertices
                     const auto fIdx = intersection.indexInInside();
-                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
@@ -457,11 +452,11 @@ public:
 
                     // find the mapped periodic vertex of all vertices on periodic boundaries
                     const auto fIdx = intersection.indexInInside();
-                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                     const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         const auto vPos = elementGeometry.corner(vIdx);
 
@@ -473,10 +468,10 @@ public:
                             if (isOutside.boundary() && isOutside.neighbor())
                             {
                                 const auto fIdxOutside = isOutside.indexInInside();
-                                const auto numFaceVertsOutside = referenceElement.size(fIdxOutside, 1, dim);
+                                const auto numFaceVertsOutside = refElement.size(fIdxOutside, 1, dim);
                                 for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
                                 {
-                                    const auto vIdxOutside = referenceElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
+                                    const auto vIdxOutside = refElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
                                     const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
                                     const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
                                     if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
-- 
GitLab


From a826c0cbe4cf6237503365e2728815da956ee408 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:43:32 +0200
Subject: [PATCH 05/10] [discretization][mpfa] replace all
 ReferenceElements::general calls with referenceElement

---
 .../discretization/cellcentered/mpfa/fvelementgeometry.hh | 6 ++----
 dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh  | 8 ++------
 dumux/discretization/cellcentered/mpfa/helper.hh          | 4 +---
 3 files changed, 5 insertions(+), 13 deletions(-)

diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index e673fff01e..9391a31ba2 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -28,7 +28,6 @@
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/iteratorrange.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/indextraits.hh>
@@ -182,7 +181,6 @@ class CCMpfaFVElementGeometry<GG, false>
     static const int dim = GridView::dimension;
     static const int dimWorld = GridView::dimensionworld;
     using CoordScalar = typename GridView::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
 public:
     //! export type of the element
@@ -377,7 +375,7 @@ private:
             const auto& e = useNeighbor ? is.outside() : element;
             const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
             const auto eg = e.geometry();
-            const auto refElement = ReferenceElements::general(eg.type());
+            const auto refElement = referenceElement(eg);
 
             // Set up a container with all relevant positions for scvf corner computation
             const auto numCorners = is.geometry().corners();
@@ -459,7 +457,7 @@ private:
             const auto& e = useNeighbor ? is.outside() : element;
             const auto indexInElement = useNeighbor ? is.indexInOutside() : is.indexInInside();
             const auto eg = e.geometry();
-            const auto refElement = ReferenceElements::general(eg.type());
+            const auto refElement = referenceElement(eg);
 
             // Set up a container with all relevant positions for scvf corner computation
             const auto numCorners = is.geometry().corners();
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 772ac838c8..4a7475dd72 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -26,8 +26,6 @@
 #ifndef DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
 #define DUMUX_DISCRETIZATION_CC_MPFA_FV_GRID_GEOMETRY_HH
 
-#include <dune/geometry/referenceelements.hh>
-
 #include <dumux/common/parameters.hh>
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/method.hh>
@@ -81,7 +79,6 @@ class CCMpfaFVGridGeometry<GV, Traits, true>
     using Intersection = typename GV::Intersection;
     using GridIndexType = typename IndexTraits<GV>::GridIndex;
     using CoordScalar = typename GV::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
     using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
 
@@ -249,7 +246,7 @@ public:
                 const auto& e = useNeighbor ? is.outside() : element;
                 const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside;
                 const auto eg = e.geometry();
-                const auto refElement = ReferenceElements::general(eg.type());
+                const auto refElement = referenceElement(eg);
 
                 // Set up a container with all relevant positions for scvf corner computation
                 const auto numCorners = is.geometry().corners();
@@ -449,7 +446,6 @@ class CCMpfaFVGridGeometry<GV, Traits, false>
     using Intersection = typename GV::Intersection;
     using GridIndexType = typename IndexTraits<GV>::GridIndex;
     using CoordScalar = typename GV::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
     using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
 
@@ -624,7 +620,7 @@ public:
                 const auto& e = useNeighbor ? is.outside() : element;
                 const auto indexInElement = useNeighbor ? is.indexInOutside() : indexInInside;
                 const auto eg = e.geometry();
-                const auto refElement = ReferenceElements::general(eg.type());
+                const auto refElement = referenceElement(eg);
 
                 // evaluate if vertices on this intersection use primary/secondary IVs
                 const bool isBranchingPoint = dim < dimWorld ? outsideIndices[indexInInside].size() > 1 : false;
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index bf04356571..f1b78cd982 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -29,7 +29,6 @@
 #include <dune/common/fmatrix.hh>
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/geometry/type.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/math.hh>
 
@@ -519,7 +518,6 @@ class CCMpfaHelper : public MpfaDimensionHelper<GridGeometry,
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
     using CoordScalar = typename GridView::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
 public:
     /*!
@@ -569,7 +567,7 @@ public:
             {
                 if (!is.neighbor() && !is.boundary())
                 {
-                    const auto refElement = ReferenceElements::general(element.type());
+                    const auto refElement = referenceElement(element);
 
                     for (int isVertex = 0; isVertex < is.geometry().corners(); ++isVertex)
                     {
-- 
GitLab


From 39a16736d5506ff451f61db34389bc24fa72017f Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:44:32 +0200
Subject: [PATCH 06/10] [multidomain][facet] replace all
 ReferenceElements::general calls with referenceElement

---
 dumux/multidomain/facet/box/couplingmapper.hh    |  4 +---
 dumux/multidomain/facet/box/fvelementgeometry.hh | 11 ++++-------
 dumux/multidomain/facet/box/fvgridgeometry.hh    | 16 ++++++----------
 dumux/multidomain/facet/codimonegridadapter.hh   |  4 +---
 dumux/multidomain/facet/enrichmenthelper.hh      |  6 ++----
 dumux/multidomain/facet/vertexmapper.hh          |  8 ++------
 6 files changed, 16 insertions(+), 33 deletions(-)

diff --git a/dumux/multidomain/facet/box/couplingmapper.hh b/dumux/multidomain/facet/box/couplingmapper.hh
index 45b2fb8a79..b54a1c5ddc 100644
--- a/dumux/multidomain/facet/box/couplingmapper.hh
+++ b/dumux/multidomain/facet/box/couplingmapper.hh
@@ -23,7 +23,6 @@
 #define DUMUX_BOX_FACETCOUPLING_MAPPER_HH
 
 #include <dune/common/indices.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/method.hh>
@@ -102,7 +101,6 @@ public:
         {
             using LowDimIndexType = typename IndexTraits<LowDimGridView>::GridIndex;
             using BulkIndexType = typename IndexTraits<BulkGridView>::GridIndex;
-            using BulkReferenceElements = Dune::ReferenceElements<typename BulkGridView::ctype, bulkDim>;
 
             const auto lowDimElemIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement);
             auto& lowDimData = this->couplingMap_(facetGridId, bulkGridId)[lowDimElemIdx];
@@ -120,7 +118,7 @@ public:
             for (auto bulkElemIdx : adjoinedEntityIndices)
             {
                 const auto bulkElement = bulkFvGridGeometry.element(bulkElemIdx);
-                const auto bulkRefElem = BulkReferenceElements::general(bulkElement.type());
+                const auto bulkRefElem = referenceElement(bulkElement);
 
                 // find the bulk element facet that lies on this low dim element (assumes conformity!)
                 bool found = false;
diff --git a/dumux/multidomain/facet/box/fvelementgeometry.hh b/dumux/multidomain/facet/box/fvelementgeometry.hh
index 910a1df660..a9477d78f8 100644
--- a/dumux/multidomain/facet/box/fvelementgeometry.hh
+++ b/dumux/multidomain/facet/box/fvelementgeometry.hh
@@ -27,7 +27,6 @@
 #include <algorithm>
 
 #include <dune/geometry/type.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/indextraits.hh>
 #include <dumux/discretization/scvandscvfiterators.hh>
@@ -58,7 +57,6 @@ class BoxFacetCouplingFVElementGeometry<GG, true>
     using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 public:
     //! export type of the element
     using Element = typename GridView::template Codim<0>::Entity;
@@ -163,7 +161,6 @@ class BoxFacetCouplingFVElementGeometry<GG, false>
 
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 
     using GeometryHelper = BoxGeometryHelper<GridView, dim,
                                              typename GG::SubControlVolume,
@@ -259,7 +256,7 @@ private:
         // get the element geometry
         auto elementGeometry = element.geometry();
         elemGeometryType_ = elementGeometry.type();
-        const auto referenceElement = ReferenceElements::general(elemGeometryType_);
+        const auto refElement = referenceElement(elementGeometry);
 
         // get the sub control volume geometries of this element
         GeometryHelper geometryHelper(elementGeometry);
@@ -283,8 +280,8 @@ private:
         for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
         {
             // find the local scv indices this scvf is connected to
-            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                         static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                         static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
             // create the sub-control volume face
             scvfs_.emplace_back(geometryHelper,
@@ -312,7 +309,7 @@ private:
 
             std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
             for (int i = 0; i < numFaceCorners; ++i)
-                vIndicesLocal[i] = static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, i, dim));
+                vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
 
             // if all vertices are living on the facet grid, this is an interiour boundary
             const bool isOnFacet = gridGeometry().isOnInteriorBoundary(element, intersection);
diff --git a/dumux/multidomain/facet/box/fvgridgeometry.hh b/dumux/multidomain/facet/box/fvgridgeometry.hh
index 04d73d6e5c..acf6463413 100644
--- a/dumux/multidomain/facet/box/fvgridgeometry.hh
+++ b/dumux/multidomain/facet/box/fvgridgeometry.hh
@@ -30,7 +30,6 @@
 #include <algorithm>
 
 #include <dune/grid/common/mcmgmapper.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/common/indextraits.hh>
@@ -104,7 +103,6 @@ class BoxFacetCouplingFVGridGeometry<Scalar, GV, true, Traits>
     static const int dim = GV::dimension;
     static const int dimWorld = GV::dimensionworld;
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
     using GeometryHelper = BoxGeometryHelper<GV, dim, typename Traits::SubControlVolume, typename Traits::SubControlVolumeFace>;
 
 public:
@@ -195,7 +193,7 @@ public:
 
             // get the element geometry
             auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // instantiate the geometry helper
             GeometryHelper geometryHelper(elementGeometry);
@@ -216,8 +214,8 @@ public:
             for (; scvfLocalIdx < element.subEntities(dim-1); ++scvfLocalIdx)
             {
                 // find the global and local scv indices this scvf is belonging to
-                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                             static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                             static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
                 // create the sub-control volume face
                 scvfs_[eIdx].emplace_back(geometryHelper,
@@ -245,7 +243,7 @@ public:
 
                 std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
                 for (int i = 0; i < numFaceCorners; ++i)
-                    vIndicesLocal[i] = static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, i, dim));
+                    vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
 
                 std::vector<LocalIndexType> gridVertexIndices(numFaceCorners);
                 for (int i = 0; i < numFaceCorners; ++i)
@@ -364,8 +362,6 @@ class BoxFacetCouplingFVGridGeometry<Scalar, GV, false, Traits>
     using Intersection = typename GV::Intersection;
     using CoordScalar = typename GV::ctype;
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-
 public:
     //! export discretization method
     static constexpr DiscretizationMethod discMethod = DiscretizationMethod::box;
@@ -449,7 +445,7 @@ public:
             numScvf_ += element.subEntities(dim-1);
 
             const auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // store the sub control volume face indices on the domain/interior boundary
             // skip handled facets (necessary for e.g. Dune::FoamGrid)
@@ -469,7 +465,7 @@ public:
 
                 std::vector<LocalIndexType> vIndicesLocal(numFaceCorners);
                 for (int i = 0; i < numFaceCorners; ++i)
-                    vIndicesLocal[i] = static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, i, dim));
+                    vIndicesLocal[i] = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, i, dim));
 
                 std::vector<GridIndexType> gridVertexIndices(numFaceCorners);
                 for (int i = 0; i < numFaceCorners; ++i)
diff --git a/dumux/multidomain/facet/codimonegridadapter.hh b/dumux/multidomain/facet/codimonegridadapter.hh
index 98be2a3be5..559391d205 100644
--- a/dumux/multidomain/facet/codimonegridadapter.hh
+++ b/dumux/multidomain/facet/codimonegridadapter.hh
@@ -30,7 +30,6 @@
 #include <memory>
 
 #include <dune/grid/common/mcmgmapper.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/indextraits.hh>
 
@@ -60,7 +59,6 @@ class CodimOneGridAdapter
     // Extract some types of the bulk grid
     using BulkGridView = typename Embeddings::template GridView<bulkGridId>;
     using BulkMapper = Dune::MultipleCodimMultipleGeomTypeMapper<BulkGridView>;
-    using BulkReferenceElements = typename Dune::ReferenceElements<typename BulkGridView::ctype, BulkGridView::dimension>;
     using BulkGridElement = typename BulkGridView::template Codim<0>::Entity;
     using BulkGridIntersection = typename BulkGridView::Intersection;
     using BulkGridVertex = typename BulkGridView::template Codim<BulkGridView::dimension>::Entity;
@@ -171,7 +169,7 @@ public:
     bool isOnFacetGrid(const BulkGridElement& element, const BulkGridIntersection& intersection) const
     {
         // Intersection lies on facet grid, if the corners of the intersection make up a facet element
-        const auto refElement = BulkReferenceElements::general(element.type());
+        const auto refElement = referenceElement(element);
         const auto numCorners = intersection.geometry().corners();
         const auto facetIdx = intersection.indexInInside();
 
diff --git a/dumux/multidomain/facet/enrichmenthelper.hh b/dumux/multidomain/facet/enrichmenthelper.hh
index 17fefcf892..40446fba71 100644
--- a/dumux/multidomain/facet/enrichmenthelper.hh
+++ b/dumux/multidomain/facet/enrichmenthelper.hh
@@ -33,7 +33,6 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/reservedvector.hh>
 
-#include <dune/geometry/referenceelements.hh>
 #include <dune/grid/common/mcmgmapper.hh>
 
 #include <dumux/common/math.hh>
@@ -61,7 +60,6 @@ class VertexEnrichmentHelper
     static_assert(dimWorld == int(CodimOneGridView::dimensionworld), "world dimension mismatch");
 
     using Intersection = typename GridView::Intersection;
-    using ReferenceElements = typename Dune::ReferenceElements<typename GridView::ctype, dim>;
     using MCMGMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
     using GridIndexType = typename IndexTraits<GridView>::GridIndex;
     using Element = typename GridView::template Codim<0>::Entity;
@@ -103,7 +101,7 @@ public:
             const auto eIdx = elementMapper.index(e);
 
             std::vector<unsigned int> handledFacets;
-            const auto refElement = ReferenceElements::general(e.type());
+            const auto refElement = referenceElement(e);
             for (const auto& is : intersections(gridView, e))
             {
                 // do not start path search on boundary intersections
@@ -273,7 +271,7 @@ private:
                 if (std::find(path.begin(), path.end(), outsideElemIdx) == path.end())
                 {
                     const auto idxInOutside = is.indexInOutside();
-                    const auto outsideRefElement = ReferenceElements::general(outsideElement.type());
+                    const auto outsideRefElement = referenceElement(outsideElement);
                     path.push_back(outsideElemIdx);
 
                     // on 2d grids, there is only going to be one more
diff --git a/dumux/multidomain/facet/vertexmapper.hh b/dumux/multidomain/facet/vertexmapper.hh
index 5b67f39076..1ad999b589 100644
--- a/dumux/multidomain/facet/vertexmapper.hh
+++ b/dumux/multidomain/facet/vertexmapper.hh
@@ -29,9 +29,7 @@
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/timer.hh>
-
 #include <dune/grid/common/mcmgmapper.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include "enrichmenthelper.hh"
 
@@ -80,10 +78,9 @@ public:
 
         // first find all bulk grid vertices on the boundary
         std::vector<bool> isOnBoundary(gridView.size(dim), false);
-        using ReferenceElements = typename Dune::ReferenceElements<typename GridView::ctype, dim>;
         for (const auto& e : elements(gridView))
         {
-            const auto refElem = ReferenceElements::general(e.type());
+            const auto refElem = referenceElement(e);
             for (const auto& is : intersections(gridView, e))
                 if (is.boundary())
                     for (int i = 0; i < is.geometry().corners(); ++i)
@@ -97,7 +94,6 @@ public:
             vertexMarkers[codimOneGridAdapter.bulkGridVertexIndex(vertex)] = true;
 
         // unmark where necessary
-        using CodimOneReferenceElements = typename Dune::ReferenceElements<typename CodimOneGridView::ctype, dim-1>;
         for (const auto& codimOneElement : elements(codimOneGridView))
         {
             // if a codimension one element has less than two embedments  we do not need to enrich
@@ -108,7 +104,7 @@ public:
             // otherwise, we check for immersed boundaries where we also must not enrich
             else
             {
-                const auto refElem = CodimOneReferenceElements::general(codimOneElement.type());
+                const auto refElem = referenceElement(codimOneElement);
                 for (const auto& intersection : intersections(codimOneGridView, codimOneElement))
                 {
                     // skip if intersection is not on boundary
-- 
GitLab


From 1931ae56e7301677edb638314f4d4a9ed3f4b583 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:44:49 +0200
Subject: [PATCH 07/10] [boxdfm] replace all ReferenceElements::general calls
 with referenceElement

---
 .../boxdfm/fvelementgeometry.hh               | 25 +++++------
 .../porousmediumflow/boxdfm/fvgridgeometry.hh | 41 ++++++++-----------
 .../porousmediumflow/boxdfm/geometryhelper.hh | 21 ++++------
 .../boxdfm/vtkoutputmodule.hh                 |  6 +--
 4 files changed, 38 insertions(+), 55 deletions(-)

diff --git a/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh b/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
index f2fd4cc350..ca0d41b767 100644
--- a/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
+++ b/dumux/porousmediumflow/boxdfm/fvelementgeometry.hh
@@ -31,7 +31,6 @@
 
 #include <dune/common/version.hh>
 #include <dune/geometry/type.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
 #include <dumux/discretization/scvandscvfiterators.hh>
@@ -62,7 +61,6 @@ class BoxDfmFVElementGeometry<GG, true>
     using Element = typename GridView::template Codim<0>::Entity;
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
 public:
     //! Export type of subcontrol volume
     using SubControlVolume = typename GG::SubControlVolume;
@@ -173,9 +171,6 @@ class BoxDfmFVElementGeometry<GG, false>
 
     using CoordScalar = typename GridView::ctype;
     using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-    using FaceReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim-1>;
-
     using GeometryHelper = BoxDfmGeometryHelper<GridView, dim,
                                                 typename GG::SubControlVolume,
                                                 typename GG::SubControlVolumeFace>;
@@ -276,7 +271,7 @@ private:
         // get the element geometry
         auto elementGeometry = element.geometry();
         elemGeometryType_ = elementGeometry.type();
-        const auto referenceElement = ReferenceElements::general(elemGeometryType_);
+        const auto refElement = referenceElement(element);
 
         // get the sub control volume geometries of this element
         GeometryHelper geometryHelper(elementGeometry);
@@ -304,8 +299,8 @@ private:
         for (; scvfLocalIdx < numInnerScvf; ++scvfLocalIdx)
         {
             // find the local scv indices this scvf is connected to
-            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                         static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                         static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
             scvfs_[scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
                                                         element,
@@ -336,7 +331,7 @@ private:
             std::vector<GridIndexType> isVertexIndices(numCorners);
             for (unsigned int vIdxLocal = 0; vIdxLocal < numCorners; ++vIdxLocal)
                 isVertexIndices[vIdxLocal] = gridGeometry().vertexMapper().subIndex(element,
-                                                                                      referenceElement.subEntity(idxInInside, 1, vIdxLocal, dim),
+                                                                                      refElement.subEntity(idxInInside, 1, vIdxLocal, dim),
                                                                                       dim);
 
             if (intersection.boundary())
@@ -344,7 +339,7 @@ private:
                 for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < numCorners; ++isScvfLocalIdx)
                 {
                     // find the scv this scvf is connected to
-                    const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, isScvfLocalIdx, dim));
+                    const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, isScvfLocalIdx, dim));
                     std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
 
                     scvfs_.emplace_back(geometryHelper,
@@ -370,7 +365,7 @@ private:
                                        intersection,
                                        isGeometry,
                                        vIdxLocal,
-                                       static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, vIdxLocal, dim)),
+                                       static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, vIdxLocal, dim)),
                                        scvLocalIdx++,
                                        idxInInside,
                                        eIdx,
@@ -379,12 +374,12 @@ private:
                 // add fracture scvf for each edge of the intersection in 3d
                 if (dim == 3)
                 {
-                    const auto& faceReferenceElement = FaceReferenceElements::general(isGeometry.type());
-                    for (unsigned int edgeIdx = 0; edgeIdx < faceReferenceElement.size(1); ++edgeIdx)
+                    const auto& faceRefElement = referenceElement(isGeometry);
+                    for (unsigned int edgeIdx = 0; edgeIdx < faceRefElement.size(1); ++edgeIdx)
                     {
                         // inside/outside scv indices in face local node numbering
-                        std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(faceReferenceElement.subEntity(edgeIdx, 1, 0, dim-1)),
-                                                                     static_cast<LocalIndexType>(faceReferenceElement.subEntity(edgeIdx, 1, 1, dim-1))});
+                        std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 0, dim-1)),
+                                                                     static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 1, dim-1))});
 
                         // add offset to get the right scv indices
                         std::for_each( localScvIndices.begin(),
diff --git a/dumux/porousmediumflow/boxdfm/fvgridgeometry.hh b/dumux/porousmediumflow/boxdfm/fvgridgeometry.hh
index 6ebb2f5874..ae89f1555e 100644
--- a/dumux/porousmediumflow/boxdfm/fvgridgeometry.hh
+++ b/dumux/porousmediumflow/boxdfm/fvgridgeometry.hh
@@ -31,7 +31,6 @@
 
 #include <unordered_map>
 
-#include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 #include <dune/geometry/multilineargeometry.hh>
 #include <dune/grid/common/mcmgmapper.hh>
@@ -108,9 +107,6 @@ class BoxDfmFVGridGeometry<Scalar, GV, true, Traits>
     static const int dimWorld = GV::dimensionworld;
     static_assert(dim == 2 || dim == 3, "The box-dfm GridGeometry is only implemented in 2 or 3 dimensions.");
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-    using FaceReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim-1>;
-
     using GeometryHelper = BoxDfmGeometryHelper<GV, dim,
                                                 typename Traits::SubControlVolume,
                                                 typename Traits::SubControlVolumeFace>;
@@ -189,7 +185,7 @@ public:
 
             // get the element geometry
             auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // instantiate the geometry helper
             GeometryHelper geometryHelper(elementGeometry);
@@ -213,8 +209,8 @@ public:
             for (; scvfLocalIdx < element.subEntities(dim-1); ++scvfLocalIdx)
             {
                 // find the global and local scv indices this scvf is belonging to
-                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
-                                                             static_cast<LocalIndexType>(referenceElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
+                std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 0, dim)),
+                                                             static_cast<LocalIndexType>(refElement.subEntity(scvfLocalIdx, dim-1, 1, dim))});
 
                 scvfs_[eIdx][scvfLocalIdx] = SubControlVolumeFace(geometryHelper,
                                                                   element,
@@ -245,7 +241,7 @@ public:
                 std::vector<GridIndexType> isVertexIndices(numCorners);
                 for (unsigned int vIdxLocal = 0; vIdxLocal < numCorners; ++vIdxLocal)
                     isVertexIndices[vIdxLocal] = this->vertexMapper().subIndex(element,
-                                                                               referenceElement.subEntity(idxInInside, 1, vIdxLocal, dim),
+                                                                               refElement.subEntity(idxInInside, 1, vIdxLocal, dim),
                                                                                dim);
                 // maybe add boundary scvf
                 if (intersection.boundary() && !intersection.neighbor())
@@ -256,7 +252,7 @@ public:
                     for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < numCorners; ++isScvfLocalIdx)
                     {
                         // find the scvs this scvf is belonging to
-                        const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, isScvfLocalIdx, dim));
+                        const LocalIndexType insideScvIdx = static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, isScvfLocalIdx, dim));
                         std::vector<LocalIndexType> localScvIndices = {insideScvIdx, insideScvIdx};
                         scvfs_[eIdx].emplace_back(geometryHelper,
                                                   intersection,
@@ -268,10 +264,10 @@ public:
 
                     // add all vertices on the intersection to the set of
                     // boundary vertices
-                    const auto numFaceVerts = referenceElement.size(idxInInside, 1, dim);
+                    const auto numFaceVerts = refElement.size(idxInInside, 1, dim);
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(idxInInside, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(idxInInside, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
@@ -295,7 +291,7 @@ public:
                                                  intersection,
                                                  isGeometry,
                                                  vIdxLocal,
-                                                 static_cast<LocalIndexType>(referenceElement.subEntity(idxInInside, 1, vIdxLocal, dim)),
+                                                 static_cast<LocalIndexType>(refElement.subEntity(idxInInside, 1, vIdxLocal, dim)),
                                                  scvLocalIdx++,
                                                  idxInInside,
                                                  eIdx,
@@ -304,12 +300,12 @@ public:
                     // add fracture scvf for each edge of the intersection in 3d
                     if (dim == 3)
                     {
-                        const auto& faceReferenceElement = FaceReferenceElements::general(isGeometry.type());
-                        for (unsigned int edgeIdx = 0; edgeIdx < faceReferenceElement.size(1); ++edgeIdx)
+                        const auto& faceRefElement = referenceElement(isGeometry);
+                        for (unsigned int edgeIdx = 0; edgeIdx < faceRefElement.size(1); ++edgeIdx)
                         {
                             // inside/outside scv indices in face local node numbering
-                            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(faceReferenceElement.subEntity(edgeIdx, 1, 0, dim-1)),
-                                                                         static_cast<LocalIndexType>(faceReferenceElement.subEntity(edgeIdx, 1, 1, dim-1))});
+                            std::vector<LocalIndexType> localScvIndices({static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 0, dim-1)),
+                                                                         static_cast<LocalIndexType>(faceRefElement.subEntity(edgeIdx, 1, 1, dim-1))});
 
                             // add offset to get the right scv indices
                             std::for_each( localScvIndices.begin(),
@@ -413,9 +409,6 @@ class BoxDfmFVGridGeometry<Scalar, GV, false, Traits>
     using Intersection = typename GV::Intersection;
     using CoordScalar = typename GV::ctype;
 
-    using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
-    using FaceReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim-1>;
-
 public:
     //! export discretization method
     static constexpr DiscretizationMethod discMethod = DiscretizationMethod::box;
@@ -482,7 +475,7 @@ public:
             numScvf_ += element.subEntities(dim-1);
 
             const auto elementGeometry = element.geometry();
-            const auto referenceElement = ReferenceElements::general(elementGeometry.type());
+            const auto refElement = referenceElement(elementGeometry);
 
             // store the sub control volume face indices on the domain boundary
             for (const auto& intersection : intersections(this->gridView(), element))
@@ -495,7 +488,7 @@ public:
                 std::vector<GridIndexType> isVertexIndices(numCorners);
                 for (unsigned int vIdxLocal = 0; vIdxLocal < numCorners; ++vIdxLocal)
                     isVertexIndices[vIdxLocal] = this->vertexMapper().subIndex(element,
-                                                                               referenceElement.subEntity(idxInInside, 1, vIdxLocal, dim),
+                                                                               refElement.subEntity(idxInInside, 1, vIdxLocal, dim),
                                                                                dim);
 
                 if (intersection.boundary() && !intersection.neighbor())
@@ -506,10 +499,10 @@ public:
                     // add all vertices on the intersection to the set of
                     // boundary vertices
                     const auto fIdx = intersection.indexInInside();
-                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto numFaceVerts = refElement.size(fIdx, 1, dim);
                     for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
                     {
-                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdx = refElement.subEntity(fIdx, 1, localVIdx, dim);
                         const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
@@ -528,7 +521,7 @@ public:
 
                     const auto isGeometry = intersection.geometry();
                     numScv_ += isGeometry.corners();
-                    numScvf_ += dim == 3 ? FaceReferenceElements::general(isGeometry.type()).size(1) : 1;
+                    numScvf_ += dim == 3 ? referenceElement(isGeometry).size(1) : 1;
                 }
             }
         }
diff --git a/dumux/porousmediumflow/boxdfm/geometryhelper.hh b/dumux/porousmediumflow/boxdfm/geometryhelper.hh
index b22242e51a..f30cdb4cb7 100644
--- a/dumux/porousmediumflow/boxdfm/geometryhelper.hh
+++ b/dumux/porousmediumflow/boxdfm/geometryhelper.hh
@@ -45,7 +45,6 @@ class BoxDfmGeometryHelper<GridView, 2, ScvType, ScvfType> : public BoxGeometryH
 
     static constexpr auto dim = GridView::dimension;
     using Scalar = typename GridView::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
 
 public:
 
@@ -67,9 +66,9 @@ public:
                    const Intersection& is,
                    unsigned int edgeIndexInIntersection = 0) const
     {
-        const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
-        const auto vIdxLocal0 = referenceElement.subEntity(is.indexInInside(), 1, 0, dim);
-        const auto vIdxLocal1 = referenceElement.subEntity(is.indexInInside(), 1, 1, dim);
+        const auto refElement = referenceElement(this->elementGeometry_);
+        const auto vIdxLocal0 = refElement.subEntity(is.indexInInside(), 1, 0, dim);
+        const auto vIdxLocal1 = refElement.subEntity(is.indexInInside(), 1, 1, dim);
         auto n = this->elementGeometry_.corner(vIdxLocal1) - this->elementGeometry_.corner(vIdxLocal0);
         n /= n.two_norm();
         return n;
@@ -88,8 +87,6 @@ class BoxDfmGeometryHelper<GridView, 3, ScvType, ScvfType> : public BoxGeometryH
     static constexpr auto dim = GridView::dimension;
     static constexpr auto dimWorld = GridView::dimensionworld;
     using Scalar = typename GridView::ctype;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
-    using FaceReferenceElements = typename Dune::ReferenceElements<Scalar, dim-1>;
 
 public:
 
@@ -101,8 +98,8 @@ public:
                                              const typename Intersection::Geometry& isGeom,
                                              unsigned int edgeIndexInIntersection) const
     {
-        const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
-        const auto faceRefElem = FaceReferenceElements::general(isGeom.type());
+        const auto refElement = referenceElement(this->elementGeometry_);
+        const auto faceRefElem = referenceElement(isGeom);
 
         // create point vector for this geometry
         typename ScvfType::Traits::GlobalPosition pi[9];
@@ -114,7 +111,7 @@ public:
         const auto idxInInside = is.indexInInside();
         for (int i = 0; i < faceRefElem.size(1); ++i)
         {
-            const auto edgeIdxLocal = referenceElement.subEntity(idxInInside, 1, i, dim-1);
+            const auto edgeIdxLocal = refElement.subEntity(idxInInside, 1, i, dim-1);
             pi[i+1] = this->p_[edgeIdxLocal+this->corners_+1];
         }
 
@@ -164,15 +161,15 @@ public:
                    const Intersection& is,
                    unsigned int edgeIndexInIntersection) const
     {
-        const auto referenceElement = ReferenceElements::general(this->elementGeometry_.type());
+        const auto refElement = referenceElement(this->elementGeometry_);
 
         // first get the intersection corners (maximum "4" is for quadrilateral face)
         typename ScvfType::Traits::GlobalPosition c[4];
 
-        const auto corners = referenceElement.size(is.indexInInside(), 1, dim);
+        const auto corners = refElement.size(is.indexInInside(), 1, dim);
         for (int i = 0; i < corners; ++i)
         {
-            const auto vIdxLocal = referenceElement.subEntity(is.indexInInside(), 1, i, dim);
+            const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, i, dim);
             c[i] = this->elementGeometry_.corner(vIdxLocal);
         }
 
diff --git a/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh b/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh
index 1b90c61a5c..a4faed52a3 100644
--- a/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh
+++ b/dumux/porousmediumflow/boxdfm/vtkoutputmodule.hh
@@ -28,7 +28,6 @@
 #include <set>
 
 #include <dune/grid/common/gridfactory.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/grid/common/mcmgmapper.hh>
 
 #include <dumux/io/vtkoutputmodule.hh>
@@ -74,7 +73,6 @@ class BoxDfmVtkOutputModule : public VtkOutputModule<GridVariables, SolutionVect
     using GridIndexType = typename GridView::IndexSet::IndexType;
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    using ReferenceElements = typename Dune::ReferenceElements<typename GridView::ctype, dim>;
 
     using Field = Vtk::template Field<GridView>;
     using FractureField = Vtk::template Field<FractureGridView>;
@@ -519,7 +517,7 @@ private:
         for (const auto& element : elements(gridView))
         {
             const auto eIdxGlobal = gridGeometry.elementMapper().index(element);
-            const auto referenceElement = ReferenceElements::general(element.type());
+            const auto refElement = referenceElement(element);
 
             for (const auto& is : intersections(gridView, element))
             {
@@ -531,7 +529,7 @@ private:
                 std::vector<GridIndexType> isVertexIndices(numCorners);
                 for (unsigned int i = 0; i < numCorners; ++i)
                     isVertexIndices[i] = gridGeometry.vertexMapper().subIndex(element,
-                                                                                referenceElement.subEntity(indexInInside, 1, i, dim),
+                                                                                refElement.subEntity(indexInInside, 1, i, dim),
                                                                                 dim);
 
                 // determine if this is a fracture facet & if it has to be inserted
-- 
GitLab


From 05ad86b66f3a83e8e111cf8d85a6515952595b26 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:45:40 +0200
Subject: [PATCH 08/10] [geometryintersection] replace all
 ReferenceElements::general calls with referenceElement

---
 dumux/common/geometry/geometryintersection.hh | 37 +++++++------------
 1 file changed, 14 insertions(+), 23 deletions(-)

diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index 306bead65f..18f611aa7f 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -27,7 +27,6 @@
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/promotiontraits.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/geometry/multilineargeometry.hh>
 
 #include <dumux/common/math.hh>
@@ -406,7 +405,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
 
 public:
     /*!
@@ -552,8 +550,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
-    using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;
 
 public:
     /*!
@@ -594,22 +590,22 @@ public:
 
             if (points.size() - numPoints1 != geo2.corners())
             {
-                const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
-                const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+                const auto refElement1 = referenceElement(geo1);
+                const auto refElement2 = referenceElement(geo2);
 
                 // add intersections of edges
                 using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
                 using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
-                for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
+                for (int i = 0; i < refElement1.size(dim1-1); ++i)
                 {
-                    const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i);
+                    const auto localEdgeGeom1 = refElement1.template geometry<dim1-1>(i);
                     const auto edge1 = SegGeometry( Dune::GeometryTypes::line,
                                                     std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)),
                                                                          geo1.global(localEdgeGeom1.corner(1))} ));
 
-                    for (int j = 0; j < referenceElement2.size(dim2-1); ++j)
+                    for (int j = 0; j < refElement2.size(dim2-1); ++j)
                     {
-                        const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j);
+                        const auto localEdgeGeom2 = refElement2.template geometry<dim2-1>(j);
                         const auto edge2 = SegGeometry( Dune::GeometryTypes::line,
                                                         std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)),
                                                                              geo2.global(localEdgeGeom2.corner(1))} ));
@@ -698,7 +694,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
 
 public:
     /*!
@@ -847,8 +842,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
-    using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;
 
 public:
     /*!
@@ -887,16 +880,16 @@ public:
         using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>;
         using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
 
-        const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
-        const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+        const auto refElement1 = referenceElement(geo1);
+        const auto refElement2 = referenceElement(geo2);
 
         // intersection policy for point-like intersections (used later)
         using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
 
         // add intersection points of all polyhedron edges (codim dim-1) with the polygon
-        for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
+        for (int i = 0; i < refElement1.size(dim1-1); ++i)
         {
-            const auto localEdgeGeom = referenceElement1.template geometry<dim1-1>(i);
+            const auto localEdgeGeom = refElement1.template geometry<dim1-1>(i);
             const auto p = geo1.global(localEdgeGeom.corner(0));
             const auto q = geo1.global(localEdgeGeom.corner(1));
             const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
@@ -908,11 +901,11 @@ public:
         }
 
         // add intersection points of all polygon faces (codim 1) with the polyhedron faces
-        for (int i = 0; i < referenceElement1.size(1); ++i)
+        for (int i = 0; i < refElement1.size(1); ++i)
         {
             const auto faceGeo = [&]()
             {
-                const auto localFaceGeo = referenceElement1.template geometry<1>(i);
+                const auto localFaceGeo = refElement1.template geometry<1>(i);
                 if (localFaceGeo.corners() == 4)
                 {
                     const auto a = geo1.global(localFaceGeo.corner(0));
@@ -932,9 +925,9 @@ public:
                 }
             }();
 
-            for (int j = 0; j < referenceElement2.size(1); ++j)
+            for (int j = 0; j < refElement2.size(1); ++j)
             {
-                const auto localEdgeGeom = referenceElement2.template geometry<1>(j);
+                const auto localEdgeGeom = refElement2.template geometry<1>(j);
                 const auto p = geo2.global(localEdgeGeom.corner(0));
                 const auto q = geo2.global(localEdgeGeom.corner(1));
 
@@ -1039,7 +1032,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
 
 public:
     /*!
@@ -1313,7 +1305,6 @@ public:
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
-    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
 
 public:
     /*!
-- 
GitLab


From 3f002dba4875e25ad4b7ed9887ce9432d740c8b7 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:59:16 +0200
Subject: [PATCH 09/10] [sequential] replace all ReferenceElements::general
 calls with referenceElement

---
 .../diffusion/cellcentered/velocity.hh        |   4 +-
 .../diffusion/cellcentered/velocity.hh        |   5 +-
 .../sequential/diffusion/mimetic/pressure.hh  |   3 +-
 .../diffusion/mimetic/pressureadaptive.hh     |   3 +-
 .../diffusion/mpfa/lmethod/2dpressure.hh      |  16 +-
 .../mpfa/lmethod/2dpressureadaptive.hh        |  20 ++-
 .../mpfa/lmethod/2dpressurevelocity.hh        |   6 +-
 .../lmethod/2dpressurevelocityadaptive.hh     |   6 +-
 .../diffusion/mpfa/lmethod/2dvelocity.hh      |  12 +-
 .../lmethod/3dinteractionvolumecontainer.hh   |  21 ++-
 .../3dinteractionvolumecontaineradaptive.hh   | 152 +++++++++---------
 .../mpfa/lmethod/3dpressurevelocity.hh        |   6 +-
 .../diffusion/mpfa/lmethod/3dvelocity.hh      |   4 +-
 .../diffusion/mpfa/omethod/2dpressure.hh      |  16 +-
 .../mpfa/omethod/2dpressurevelocity.hh        |   7 +-
 .../diffusion/mpfa/omethod/2dvelocity.hh      |  13 +-
 .../2p2c/sequential/fv3dpressureadaptive.hh   |  39 +++--
 17 files changed, 151 insertions(+), 182 deletions(-)

diff --git a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/velocity.hh b/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/velocity.hh
index d5fd72e4bb..d33d269af4 100644
--- a/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/velocity.hh
+++ b/dumux/porousmediumflow/1p/sequential/diffusion/cellcentered/velocity.hh
@@ -27,7 +27,6 @@
 
 #include <dumux/porousmediumflow/1p/sequential/properties.hh>
 
-
 namespace Dumux {
 /*!
  * \ingroup SequentialOnePModel
@@ -126,8 +125,7 @@ public:
             const typename Element::Geometry& geometry = element.geometry();
 
             // get corresponding reference element
-            using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-            const auto refElement = ReferenceElements::general(geometry.type());
+            const auto refElement = referenceElement(geometry);
 
             const int numberOfFaces = refElement.size(1);
             std::vector<Scalar> flux(numberOfFaces,0);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
index 175fcea4ec..b4aa4760ce 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
@@ -87,8 +87,6 @@ class FVVelocity2P
         dim = GridView::dimension, dimWorld = GridView::dimensionworld
     };
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     enum
     {
         pw = Indices::pressureW,
@@ -203,8 +201,7 @@ public:
                 const typename Element::Geometry& geometry = element.geometry();
 
                 // get corresponding reference element
-                using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-                const auto refElement = ReferenceElements::general(geometry.type());
+                const auto refElement = referenceElement(geometry);
 
                 const int numberOfFaces=refElement.size(1);
                 std::vector<Scalar> fluxW(numberOfFaces,0);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
index 8176c2be08..a03080ad98 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
@@ -294,8 +294,7 @@ public:
                 const typename Element::Geometry& geometry = element.geometry();
 
                 // get corresponding reference element
-                using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-                const auto refElement = ReferenceElements::general(geometry.type());
+                const auto refElement = referenceElement(geometry);
 
                 const int numberOfFaces=refElement.size(1);
                 std::vector<Scalar> fluxW(numberOfFaces,0);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
index ba85227d09..e89c758494 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
@@ -308,8 +308,7 @@ public:
                 const typename Element::Geometry& geometry = element.geometry();
 
                 // get corresponding reference element
-                using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-                const auto refElement = ReferenceElements::general(geometry.type());
+                const auto refElement = referenceElement(geometry);
 
                 const int numberOfFaces=refElement.size(1);
                 std::vector<Scalar> fluxW(numberOfFaces,0);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
index d7e91e2a1d..901e085612 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
@@ -81,8 +81,6 @@ class FvMpfaL2dPressure2p: public FVPressure<TypeTag>
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -733,7 +731,7 @@ void FvMpfaL2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
 
             // get the intersection node /bar^{x_3} between 'intersection12'
             // and 'intersection14', denoted as 'corner1234'
-            const auto referenceElement = ReferenceElements::general(element.type());
+            const auto refElement = referenceElement(element);
 
             GlobalPosition corner1234(0);
 
@@ -744,13 +742,13 @@ void FvMpfaL2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
             {
                 bool finished = false;
 
-                int localVertIdx12corner = referenceElement.subEntity(indexInInside12, 1, i, dim);
+                int localVertIdx12corner = refElement.subEntity(indexInInside12, 1, i, dim);
 
                 int globalVertIdx12corner = problem_.variables().vertexMapper().subIndex(element, localVertIdx12corner, dim);
 
                 for (int j = 0; j < intersection14.geometry().corners(); ++j)
                 {
-                    int localVertIdx14corner = referenceElement.subEntity(indexInInside14, 1, j, dim);
+                    int localVertIdx14corner = refElement.subEntity(indexInInside14, 1, j, dim);
 
                     int globalVertIdx14corner = problem_.variables().vertexMapper().subIndex(element, localVertIdx14corner, dim);
 
@@ -961,7 +959,7 @@ void FvMpfaL2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection2.geometry().corners(); ++i)
                             {
-                                int localVertIdx2corner = referenceElement.subEntity(intersection2.indexInInside(), dim - 1, i,
+                                int localVertIdx2corner = refElement.subEntity(intersection2.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx2corner = problem_.variables().index(
@@ -1092,7 +1090,7 @@ void FvMpfaL2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection4.geometry().corners(); ++i)
                             {
-                                int localVertIdx4corner = referenceElement.subEntity(intersection4.indexInInside(), dim - 1, i,
+                                int localVertIdx4corner = refElement.subEntity(intersection4.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx4corner = problem_.variables().index(
@@ -1631,9 +1629,9 @@ void FvMpfaL2dPressure2p<TypeTag>::assemble()
                         {
                             int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                            const auto referenceElement = ReferenceElements::general(element.type());
+                            const auto refElement = referenceElement(element);
 
-                            const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                            const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                             const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
index 8ed095778b..137e313910 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
@@ -84,8 +84,6 @@ class FvMpfaL2dPressure2pAdaptive: public FVPressure<TypeTag>
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -776,7 +774,7 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
         // get index
         int eIdxGlobal1 = problem_.variables().index(element);
 
-        const auto referenceElement = ReferenceElements::general(element.type());
+        const auto refElement = referenceElement(element);
 
         const auto isEndIt12 = problem_.gridView().iend(element);
         for (auto isIt12 = problem_.gridView().ibegin(element); isIt12 != isEndIt12; ++isIt12)
@@ -807,13 +805,13 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
                 // get the global coordinate and global vertex index of corner1234
                 for (int i = 0; i < intersection12.geometry().corners(); ++i)
                 {
-                    int localVertIdx12corner = referenceElement.subEntity(indexInInside12, dim - 1, i, dim);
+                    int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
 
                     int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
 
                     for (int j = 0; j < intersection14.geometry().corners(); ++j)
                     {
-                        int localVertIdx14corner = referenceElement.subEntity(indexInInside14, dim - 1, j, dim);
+                        int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
 
                         int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
 
@@ -1245,7 +1243,7 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection2.geometry().corners(); ++i)
                             {
-                                int localVertIdx2corner = referenceElement.subEntity(intersection2.indexInInside(), dim - 1, i,
+                                int localVertIdx2corner = refElement.subEntity(intersection2.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx2corner = problem_.variables().index(element2.template subEntity<dim>(localVertIdx2corner));
@@ -1334,13 +1332,13 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
                 {
                     bool finished = false;
 
-                    int localVertIdx12corner = referenceElement.subEntity(indexInInside12, dim - 1, i, dim);
+                    int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
 
                     int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
 
                     for (int j = 0; j < intersection14.geometry().corners(); ++j)
                     {
-                        int localVertIdx14corner = referenceElement.subEntity(indexInInside14, dim - 1, j, dim);
+                        int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
 
                         int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
 
@@ -1478,7 +1476,7 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection4.geometry().corners(); ++i)
                             {
-                                int localVertIdx4corner = referenceElement.subEntity(intersection4.indexInInside(), dim - 1, i,
+                                int localVertIdx4corner = refElement.subEntity(intersection4.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx4corner = problem_.variables().index(element4.template subEntity<dim>(localVertIdx4corner));
@@ -2350,9 +2348,9 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::assemble()
                         {
                             int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                            const auto referenceElement = ReferenceElements::general(element.type());
+                            const auto refElement = referenceElement(element);
 
-                            const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                            const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                             const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
index a3f9de827d..6047876798 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
@@ -60,8 +60,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2p: public FvMpfaL2dPress
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -306,7 +304,7 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
 
-    const auto referenceElement = ReferenceElements::general(elementI.type());
+    const auto refElement = referenceElement(elementI);
 
     int indexInInside = intersection.indexInInside();
     int indexInOutside = intersection.indexInOutside();
@@ -315,7 +313,7 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     for (int vIdx = 0; vIdx < numVertices; vIdx++)
     {
-        int localVertIdx = referenceElement.subEntity(indexInInside, dim - 1, vIdx, dim);
+        int localVertIdx = refElement.subEntity(indexInInside, dim - 1, vIdx, dim);
 
         int vIdxGlobal = problem_.variables().index(elementI.template subEntity<dim>(localVertIdx));
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
index a4f86fdbdf..1cee737ca0 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
@@ -62,8 +62,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2pAdaptive: public FvMpfa
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -344,7 +342,7 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Inter
 
     CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
 
-    const auto referenceElement = ReferenceElements::general(elementI.type());
+    const auto refElement = referenceElement(elementI);
 
     int indexInInside = intersection.indexInInside();
     int indexInOutside = intersection.indexInOutside();
@@ -372,7 +370,7 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Inter
 
     for (int vIdx = 0; vIdx < numVertices; vIdx++)
     {
-        int localVertIdx = referenceElement.subEntity(fIdx, dim - 1, vIdx, dim);
+        int localVertIdx = refElement.subEntity(fIdx, dim - 1, vIdx, dim);
 
         int vIdxGlobal = 0;
                 if (levelI >= levelJ)
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
index a521a3b319..f53bbe6510 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
@@ -65,8 +65,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2p
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -209,7 +207,7 @@ public:
                 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
                 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
 
-                const DimVector& localPos = ReferenceElements::general(element.type()).position(0, 0);
+                const DimVector& localPos = referenceElement(element).position(0, 0);
 
                 // get the transposed Jacobian of the element mapping
                 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
@@ -755,9 +753,9 @@ void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
                 {
                     int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                    const auto referenceElement = ReferenceElements::general(element.type());
+                    const auto refElement = referenceElement(element);
 
-                    const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                    const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                     const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
@@ -864,9 +862,9 @@ void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
                 {
                     int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                    const auto referenceElement = ReferenceElements::general(element.type());
+                    const auto refElement = referenceElement(element);
 
-                    const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                    const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                     const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontainer.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontainer.hh
index fdc011c85b..882a5da5ef 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontainer.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontainer.hh
@@ -55,8 +55,6 @@ class FvMpfaL3dInteractionVolumeContainer
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
     using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
@@ -376,8 +374,7 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeIntersectionInfo(const E
     int eIdxGlobal = problem_.variables().index(element);
 
     const ElementGeometry& geometry = element.geometry();
-
-    const auto referenceElement = ReferenceElements::general(geometry.type());
+    const auto refElement = referenceElement(geometry);
 
     int levelI = element.level();
 
@@ -417,7 +414,7 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeIntersectionInfo(const E
 
             for (int i = 0; i < isGeometry.corners(); i++)
             {
-                int localVertIdx = referenceElement.subEntity(indexInInside, 1, i, dim);
+                int localVertIdx = refElement.subEntity(indexInInside, 1, i, dim);
 
                 int vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, localVertIdx, dim);
 
@@ -1353,20 +1350,20 @@ void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInnerInteractionVolume(I
         const ElementGeometry& geometry1 = element1.geometry();
         const ElementGeometry& geometry8 = element8.geometry();
 
-        const auto referenceElement = ReferenceElements::general(geometry1.type());
+        const auto refElement = referenceElement(geometry1);
 
-        DimVector edgeCoord(geometry1.global(referenceElement.position(9, dim - 1)));
+        DimVector edgeCoord(geometry1.global(refElement.position(9, dim - 1)));
         interactionVolume.setEdgePosition(edgeCoord, 2);
-        edgeCoord = geometry1.global(referenceElement.position(3, dim - 1));
+        edgeCoord = geometry1.global(refElement.position(3, dim - 1));
         interactionVolume.setEdgePosition(edgeCoord, 0);
-        edgeCoord = geometry1.global(referenceElement.position(11, dim - 1));
+        edgeCoord = geometry1.global(refElement.position(11, dim - 1));
         interactionVolume.setEdgePosition(edgeCoord, 5);
 
-        edgeCoord = geometry8.global(referenceElement.position(4, dim - 1));
+        edgeCoord = geometry8.global(refElement.position(4, dim - 1));
         interactionVolume.setEdgePosition(edgeCoord, 4);
-        edgeCoord = geometry8.global(referenceElement.position(6, dim - 1));
+        edgeCoord = geometry8.global(refElement.position(6, dim - 1));
         interactionVolume.setEdgePosition(edgeCoord, 3);
-        edgeCoord = geometry8.global(referenceElement.position(0, dim - 1));
+        edgeCoord = geometry8.global(refElement.position(0, dim - 1));
         interactionVolume.setEdgePosition(edgeCoord, 1);
     }
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontaineradaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontaineradaptive.hh
index 5e42bb79ee..aca3638b79 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontaineradaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dinteractionvolumecontaineradaptive.hh
@@ -57,8 +57,6 @@ class FvMpfaL3dInteractionVolumeContainerAdaptive: public FvMpfaL3dInteractionVo
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
     using Element = typename GridView::Traits::template Codim<0>::Entity;
@@ -188,94 +186,94 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteraction
 
             const ElementGeometry& geometry = element.geometry();
 
-            const auto referenceElement = ReferenceElements::general(geometry.type());
+            const auto refElement = referenceElement(geometry);
 
             switch (idx)
             {
             case 0:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(9, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 2);
-                    edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(3, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 0);
-                    edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(11, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 5);
 
                     break;
                 }
             case 1:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 0);
-                    edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(8, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 2);
-                    edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(11, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 3);
 
                     break;
                 }
             case 2:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 0);
-                    edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(9, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 4);
-                    edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(10, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 5);
 
                     break;
                 }
             case 3:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 0);
-                    edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(8, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 4);
-                    edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(10, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 3);
 
                     break;
                 }
             case 4:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(3, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 1);
-                    edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(5, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 2);
-                    edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(7, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 5);
 
                     break;
                 }
             case 5:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 1);
-                    edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(4, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 2);
-                    edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(7, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 3);
 
                     break;
                 }
             case 6:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 1);
-                    edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(5, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 4);
-                    edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(6, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 5);
 
                     break;
                 }
             case 7:
                 {
-                    DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
+                    DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
                     interactionVolume.setEdgePosition(edgeCoord, 1);
-                    edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(4, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 4);
-                    edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                    edgeCoord = geometry.global(refElement.position(6, dim - 1));
                     interactionVolume.setEdgePosition(edgeCoord, 3);
 
                     break;
@@ -349,94 +347,94 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
         const ElementGeometry& geometry = element.geometry();
 
-        const auto referenceElement = ReferenceElements::general(geometry.type());
+        const auto refElement = referenceElement(geometry);
 
         switch (idx)
         {
         case 0:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(9, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 2);
-                edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
+                edgeCoord = geometry.global(refElement.position(3, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 0);
-                edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                edgeCoord = geometry.global(refElement.position(11, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 5);
 
                 break;
             }
         case 1:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 0);
-                edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                edgeCoord = geometry.global(refElement.position(8, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 2);
-                edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                edgeCoord = geometry.global(refElement.position(11, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 3);
 
                 break;
             }
         case 2:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 0);
-                edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
+                edgeCoord = geometry.global(refElement.position(9, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 4);
-                edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                edgeCoord = geometry.global(refElement.position(10, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 5);
 
                 break;
             }
         case 3:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 0);
-                edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                edgeCoord = geometry.global(refElement.position(8, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 4);
-                edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                edgeCoord = geometry.global(refElement.position(10, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 3);
 
                 break;
             }
         case 4:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(3, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 1);
-                edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                edgeCoord = geometry.global(refElement.position(5, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 2);
-                edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                edgeCoord = geometry.global(refElement.position(7, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 5);
 
                 break;
             }
         case 5:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 1);
-                edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                edgeCoord = geometry.global(refElement.position(4, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 2);
-                edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                edgeCoord = geometry.global(refElement.position(7, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 3);
 
                 break;
             }
         case 6:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 1);
-                edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                edgeCoord = geometry.global(refElement.position(5, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 4);
-                edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                edgeCoord = geometry.global(refElement.position(6, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 5);
 
                 break;
             }
         case 7:
             {
-                DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
+                DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
                 interactionVolume.setEdgePosition(edgeCoord, 1);
-                edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                edgeCoord = geometry.global(refElement.position(4, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 4);
-                edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                edgeCoord = geometry.global(refElement.position(6, dim - 1));
                 interactionVolume.setEdgePosition(edgeCoord, 3);
 
                 break;
@@ -1118,7 +1116,7 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
 
                 const ElementGeometry& geometry = element5.geometry();
 
-                const auto referenceElement = ReferenceElements::general(geometry.type());
+                const auto refElement = referenceElement(geometry);
 
                 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
                 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
@@ -1131,13 +1129,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 2:
-                            edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(9, dim - 1));
                             break;
                         case 0:
-                            edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(3, dim - 1));
                             break;
                         case 5:
-                            edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(11, dim - 1));
                             break;
                         }
 
@@ -1148,13 +1146,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 0:
-                            edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(2, dim - 1));
                             break;
                         case 2:
-                            edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(8, dim - 1));
                             break;
                         case 3:
-                            edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(11, dim - 1));
                             break;
                         }
 
@@ -1165,13 +1163,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 0:
-                            edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(1, dim - 1));
                             break;
                         case 4:
-                            edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(9, dim - 1));
                             break;
                         case 5:
-                            edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(10, dim - 1));
                             break;
                         }
 
@@ -1182,13 +1180,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 0:
-                            edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(0, dim - 1));
                             break;
                         case 4:
-                            edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(8, dim - 1));
                             break;
                         case 3:
-                            edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(10, dim - 1));
                             break;
                         }
 
@@ -1199,13 +1197,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 1:
-                            edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(3, dim - 1));
                             break;
                         case 2:
-                            edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(5, dim - 1));
                             break;
                         case 5:
-                            edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(7, dim - 1));
                             break;
                         }
 
@@ -1216,13 +1214,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 1:
-                            edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(2, dim - 1));
                             break;
                         case 2:
-                            edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(4, dim - 1));
                             break;
                         case 3:
-                            edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(7, dim - 1));
                             break;
                         }
 
@@ -1233,13 +1231,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 1:
-                            edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(1, dim - 1));
                             break;
                         case 4:
-                            edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(5, dim - 1));
                             break;
                         case 5:
-                            edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(6, dim - 1));
                             break;
                         }
 
@@ -1250,13 +1248,13 @@ void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInter
                         switch (oldEdgeIdx)
                         {
                         case 1:
-                            edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(0, dim - 1));
                             break;
                         case 4:
-                            edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(4, dim - 1));
                             break;
                         case 3:
-                            edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
+                            edgeCoord = geometry.global(refElement.position(6, dim - 1));
                             break;
                         }
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
index b0db69570b..2cbdee37a4 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
@@ -79,8 +79,6 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2p: public FvMpfaL3dPress
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
     using Intersection = typename GridView::Intersection;
 
@@ -307,7 +305,7 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
 
-    const auto referenceElement = ReferenceElements::general(elementI.type());
+    const auto refElement = referenceElement(elementI);
 
     int indexInInside = intersection.indexInInside();
     int indexInOutside = intersection.indexInOutside();
@@ -316,7 +314,7 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     for (int vIdx = 0; vIdx < numVertices; vIdx++)
     {
-        int localVertIdx = referenceElement.subEntity(indexInInside, 1, vIdx, dim);
+        int localVertIdx = refElement.subEntity(indexInInside, 1, vIdx, dim);
 
         int vIdxGlobal = problem_.variables().index(elementI.template subEntity<dim>(localVertIdx));
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
index 39d7cb2b49..0eee350450 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
@@ -62,8 +62,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2p
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -222,7 +220,7 @@ public:
                 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
                 refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
 
-                const DimVector& localPos = ReferenceElements::general(element.type()).position(0, 0);
+                const DimVector& localPos = referenceElement(element).position(0, 0);
 
                 // get the transposed Jacobian of the element mapping
                 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
index 1271b3a466..c45cdb41f3 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
@@ -76,8 +76,6 @@ class FvMpfaO2dPressure2p: public FVPressure<TypeTag>
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -729,7 +727,7 @@ void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
             // get the intersection node /bar^{x_3} between 'isIt12' and 'isIt14', denoted as 'corner1234'
             // initialization of corner1234
 
-            const auto referenceElement = ReferenceElements::general(element.type());
+            const auto refElement = referenceElement(element);
 
             GlobalPosition corner1234(0);
 
@@ -742,13 +740,13 @@ void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
 
                 const GlobalPosition& isIt12corner = intersection12.geometry().corner(i);
 
-                int localVertIdx12corner = referenceElement.subEntity(indexInInside12, dim - 1, i, dim);
+                int localVertIdx12corner = refElement.subEntity(indexInInside12, dim - 1, i, dim);
 
                 int globalVertIdx12corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx12corner));
 
                 for (int j = 0; j < intersection14.geometry().corners(); ++j)
                 {
-                    int localVertIdx14corner = referenceElement.subEntity(indexInInside14, dim - 1, j, dim);
+                    int localVertIdx14corner = refElement.subEntity(indexInInside14, dim - 1, j, dim);
 
                     int globalVertIdx14corner = problem_.variables().index(element.template subEntity<dim>(localVertIdx14corner));
 
@@ -1018,7 +1016,7 @@ void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection2.geometry().corners(); ++i)
                             {
-                                int localVertIdx2corner = referenceElement.subEntity(intersection2.indexInInside(), dim - 1, i,
+                                int localVertIdx2corner = refElement.subEntity(intersection2.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx2corner = problem_.variables().index(element2.template subEntity<dim>(localVertIdx2corner));
@@ -1164,7 +1162,7 @@ void FvMpfaO2dPressure2p<TypeTag>::storeInteractionVolumeInfo()
                         {
                             for (int i = 0; i < intersection4.geometry().corners(); ++i)
                             {
-                                int localVertIdx4corner = referenceElement.subEntity(intersection4.indexInInside(), dim - 1, i,
+                                int localVertIdx4corner = refElement.subEntity(intersection4.indexInInside(), dim - 1, i,
                                         dim);
 
                                 int globalVertIdx4corner = problem_.variables().index(element4.template subEntity<dim>(localVertIdx4corner));
@@ -1730,9 +1728,9 @@ void FvMpfaO2dPressure2p<TypeTag>::assemble()
                         {
                             int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                            const auto referenceElement = ReferenceElements::general(element.type());
+                            const auto refElement = referenceElement(element);
 
-                            const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                            const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                             const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
index 410570eb85..2c2fb4a3d8 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
@@ -23,6 +23,7 @@
 #define DUMUX_MPFAO2DPRESSUREVELOCITIES2P_HH
 
 #include <dune/common/float_cmp.hh>
+
 #include "2dpressure.hh"
 #include "2dvelocity.hh"
 
@@ -67,8 +68,6 @@ template<class TypeTag> class FvMpfaO2dPressureVelocity2p: public FvMpfaO2dPress
     using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
     using Intersection = typename GridView::Intersection;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using InteractionVolume = FVMPFAOInteractionVolume<TypeTag>;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -304,7 +303,7 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
 
-    const auto referenceElement = ReferenceElements::general(elementI.type());
+    const auto refElement = referenceElement(elementI);
 
     int indexInInside = intersection.indexInInside();
     int indexInOutside = intersection.indexInOutside();
@@ -313,7 +312,7 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocity(const Intersection&
 
     for (int vIdx = 0; vIdx < numVertices; vIdx++)
     {
-        int localVertIdx = referenceElement.subEntity(indexInInside, dim - 1, vIdx, dim);
+        int localVertIdx = refElement.subEntity(indexInInside, dim - 1, vIdx, dim);
 
         int vIdxGlobal = problem_.variables().index(elementI.template subEntity<dim>(localVertIdx));
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
index d9b802d68a..25a3296cd7 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
@@ -25,6 +25,7 @@
 #define DUMUX_FVMPFAO2DVELOCITY2P_HH
 
 #include <dune/grid/common/gridenums.hh>
+
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
@@ -63,8 +64,6 @@ template<class TypeTag> class FvMpfaO2dVelocity2P
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
     using MaterialLaw = typename SpatialParams::MaterialLaw;
 
@@ -206,7 +205,7 @@ public:
                 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
                 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
 
-                const DimVector& localPos = ReferenceElements::general(element.type()).position(0, 0);
+                const DimVector& localPos = referenceElement(element).position(0, 0);
 
                 // get the transposed Jacobian of the element mapping
                 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
@@ -587,9 +586,9 @@ void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
                 {
                     int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                    const auto referenceElement = ReferenceElements::general(element.type());
+                    const auto refElement = referenceElement(element);
 
-                    const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                    const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                     const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
@@ -696,9 +695,9 @@ void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
                 {
                     int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
 
-                    const auto referenceElement = ReferenceElements::general(element.type());
+                    const auto refElement = referenceElement(element);
 
-                    const LocalPosition& localPos = referenceElement.position(boundaryFaceIdx, 1);
+                    const LocalPosition& localPos = refElement.position(boundaryFaceIdx, 1);
 
                     const GlobalPosition& globalPosFace = element.geometry().global(localPos);
 
diff --git a/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh b/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
index e4972f7d46..2b0248b33a 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
@@ -126,7 +126,6 @@ template<class TypeTag> class FV3dPressure2P2CAdaptive
     // using declarations to abbreviate several dune classes...
     using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
     using Element = typename GridView::Traits::template Codim<0>::Entity;
-    using ReferenceElementContainer = Dune::ReferenceElements<Scalar, dim>;
 
     using Grid = typename GridView::Grid;
     using Intersection = typename GridView::Intersection;
@@ -1436,21 +1435,21 @@ int FV3dPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const Intersect
     int localFace24 = face24->indexInInside();
     int localFace26 = face26->indexInInside();
 
-    const auto referenceElement = ReferenceElementContainer::general(neighbor.type());
+    const auto refElement = referenceElement(neighbor);
 
     //find 'x'5 = edgeCoord1226
     int edge1226;
     // search through edges of face 12
-    for(int nectarine=0; nectarine < referenceElement.size(localFace12, 1, dim-1); nectarine++)
+    for(int nectarine=0; nectarine < refElement.size(localFace12, 1, dim-1); nectarine++)
     {
         // get local Idx of edge on face 12
-        int localEdgeOn12 = referenceElement.subEntity(localFace12, 1, nectarine, dim-1);
+        int localEdgeOn12 = refElement.subEntity(localFace12, 1, nectarine, dim-1);
         // search through edges of face 26
-        for(int plum = 0; plum < referenceElement.size(localFace26, 1,dim-1); plum++)
+        for(int plum = 0; plum < refElement.size(localFace26, 1,dim-1); plum++)
         {
-//            int localEdge26 = referenceElement.subEntity(localFace26, 1, plum, dim-1);
-            if(referenceElement.subEntity(localFace12, 1, nectarine, dim-1)
-                    == referenceElement.subEntity(localFace26, 1, plum, dim-1))
+//            int localEdge26 = refElement.subEntity(localFace26, 1, plum, dim-1);
+            if(refElement.subEntity(localFace12, 1, nectarine, dim-1)
+                    == refElement.subEntity(localFace26, 1, plum, dim-1))
             {
                 edge1226 = localEdgeOn12;
                 break;
@@ -1458,19 +1457,19 @@ int FV3dPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const Intersect
         }
     }
     GlobalPosition edgeCoord1226 =  // 'x'5
-            neighbor.geometry().global(referenceElement.position(edge1226, dim-1));
+            neighbor.geometry().global(refElement.position(edge1226, dim-1));
 
     //find 'x'4 = edgeCoord1224
     int edge1224;
     // search through edges of face 12
-    for(int nectarine=0; nectarine < referenceElement.size(localFace12, 1, dim-1); nectarine++)
+    for(int nectarine=0; nectarine < refElement.size(localFace12, 1, dim-1); nectarine++)
     {
         // get local Idx of edge on face 12
-        int localEdgeOn12 = referenceElement.subEntity(localFace12, 1, nectarine, dim-1);
+        int localEdgeOn12 = refElement.subEntity(localFace12, 1, nectarine, dim-1);
         // search through edges of face 24
-        for(int plum = 0; plum < referenceElement.size(localFace24, 1, dim-1); plum++)
+        for(int plum = 0; plum < refElement.size(localFace24, 1, dim-1); plum++)
         {
-            int localEdge24 = referenceElement.subEntity(localFace24, 1, plum, dim-1);
+            int localEdge24 = refElement.subEntity(localFace24, 1, plum, dim-1);
             if(localEdgeOn12 == localEdge24)
             {
                 edge1224 = localEdgeOn12;
@@ -1479,19 +1478,19 @@ int FV3dPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const Intersect
         }
     }
     GlobalPosition edgeCoord1224 =  // 'x'4
-            neighbor.geometry().global(referenceElement.position(edge1224, dim-1));
+            neighbor.geometry().global(refElement.position(edge1224, dim-1));
 
     //find 'x'6 = edgeCoord2426
     int edge2426;
     // search through edges of face 24
-    for(int nectarine=0; nectarine < referenceElement.size(localFace24, 1, dim-1); nectarine++)
+    for(int nectarine=0; nectarine < refElement.size(localFace24, 1, dim-1); nectarine++)
     {
         // get local Idx of edge on face 24
-        int localEdgeOn24 = referenceElement.subEntity(localFace24, 1, nectarine, dim-1);
+        int localEdgeOn24 = refElement.subEntity(localFace24, 1, nectarine, dim-1);
         // search through edges of face 26
-        for(int plum = 0; plum < referenceElement.size(localFace26, 1, dim-1); plum++)
+        for(int plum = 0; plum < refElement.size(localFace26, 1, dim-1); plum++)
         {
-            int localEdge26 = referenceElement.subEntity(localFace26, 1, plum, dim-1);
+            int localEdge26 = refElement.subEntity(localFace26, 1, plum, dim-1);
             if(localEdgeOn24 == localEdge26)
             {
                 edge2426 = localEdgeOn24;
@@ -1500,7 +1499,7 @@ int FV3dPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const Intersect
         }
     }
     GlobalPosition edgeCoord2426 =   // 'x'6
-            neighbor.geometry().global(referenceElement.position(edge2426, dim-1));
+            neighbor.geometry().global(refElement.position(edge2426, dim-1));
 
     /** 2) Calculate omega, chi for matrices  **/
     // center of face in global coordinates, i.e., the midpoint of face 'isIt24'
@@ -1746,7 +1745,7 @@ int FV3dPressure2P2CAdaptive<TypeTag>::computeTransmissibilities(const Intersect
 
                 // get postion as seen from element
                 GlobalPosition vertexOnElement
-                    = referenceElement.position(verticeSmall, dim);
+                    = refElement.position(verticeSmall, dim);
 
                 for (int indexOnFace = 0; indexOnFace < 4; indexOnFace++)
                 {
-- 
GitLab


From da64056bd6306040eba2268bf4cf20ab07c6454c Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Thu, 28 May 2020 10:47:00 +0200
Subject: [PATCH 10/10] [test] replace all ReferenceElements::general calls
 with referenceElement

---
 test/io/gridmanager/gridmanagertests.hh       |  4 +--
 .../1p/sequential/resultevaluation.hh         |  4 +--
 .../1p/sequential/resultevaluation3d.hh       | 29 ++++---------------
 3 files changed, 7 insertions(+), 30 deletions(-)

diff --git a/test/io/gridmanager/gridmanagertests.hh b/test/io/gridmanager/gridmanagertests.hh
index 7dbbc38f24..9ba89f6aa4 100644
--- a/test/io/gridmanager/gridmanagertests.hh
+++ b/test/io/gridmanager/gridmanagertests.hh
@@ -25,7 +25,6 @@
 #define DUMUX_TEST_IO_GRIDMANAGER_TESTS_HH
 
 #include <dune/common/version.hh>
-#include <dune/geometry/referenceelements.hh>
 #include <dune/grid/common/datahandleif.hh>
 #include <dune/grid/common/mcmgmapper.hh>
 #include <dune/grid/io/file/vtk.hh>
@@ -106,7 +105,6 @@ class GridManagerTests
     using Scalar = double;
     static const int dim = Grid::dimension;
     using GridManager = typename Dumux::GridManager<Grid>;
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
 
 public:
 
@@ -261,7 +259,7 @@ private:
                 if (!intersection.boundary())
                     continue;
 
-                const auto refElement = ReferenceElements::general(element.geometry().type());
+                const auto refElement = referenceElement(element.geometry());
                 const auto numVertices = refElement.size(intersection.indexInInside(), 1, dim);
                 // loop over vertices of the intersection facet
                 for(std::size_t vIdx = 0; vIdx < numVertices; vIdx++)
diff --git a/test/porousmediumflow/1p/sequential/resultevaluation.hh b/test/porousmediumflow/1p/sequential/resultevaluation.hh
index 775c890260..5c24ed639f 100644
--- a/test/porousmediumflow/1p/sequential/resultevaluation.hh
+++ b/test/porousmediumflow/1p/sequential/resultevaluation.hh
@@ -102,9 +102,7 @@ public:
         {
             // element geometry
             const Geometry& geometry = element.geometry();
-            Dune::GeometryType geomType = geometry.type();
-            using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-            const Dune::FieldVector<Scalar,dim>& local = ReferenceElements::general(geomType).position(0, 0);
+            const Dune::FieldVector<Scalar,dim>& local = referenceElement(geometry).position(0, 0);
             Dune::FieldVector<Scalar,dim> globalPos = geometry.global(local);
 
             Scalar volume = geometry.volume();
diff --git a/test/porousmediumflow/1p/sequential/resultevaluation3d.hh b/test/porousmediumflow/1p/sequential/resultevaluation3d.hh
index 2c4b22df8a..b3131fa96f 100644
--- a/test/porousmediumflow/1p/sequential/resultevaluation3d.hh
+++ b/test/porousmediumflow/1p/sequential/resultevaluation3d.hh
@@ -92,7 +92,6 @@ public:
         enum{dim = Grid::dimension};
 
         using JacobianInverseTransposed = typename Geometry::JacobianInverseTransposed;
-        using ReferenceElements = Dune::ReferenceElements<ct, dim>;
         const GV& gridview(grid.levelGridView(grid.maxLevel()));
         const IS& indexset(gridview.indexSet());
         Mapper elementMapper(gridview, Dune::mcmgElementLayout());
@@ -159,12 +158,8 @@ public:
             // element geometry
             const Geometry& geometry = element.geometry();
 
-            // cell geometry type
-            Dune::GeometryType gt = geometry.type();
-
             // cell center in reference element
-            const Dune::FieldVector<ct,dim>&
-                local = ReferenceElements::general(gt).position(0,0);
+            const Dune::FieldVector<ct,dim>& local = referenceElement(geometry).position(0,0);
 
             // get global coordinate of cell center
             Dune::FieldVector<ct,dim> globalPos = geometry.global(local);
@@ -425,12 +420,10 @@ public:
         using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
         using SolVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
         using JacobianInverseTransposed = typename Geometry::JacobianInverseTransposed;
-        using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
 
         ElementMapper elementMapper(gridView, Dune::mcmgElementLayout());
         SolVector exactSol(gridView.size(0));
 
-
         uMinExact = 1e100;
         uMaxExact = -1e100;
         uMin = 1e100;
@@ -469,10 +462,7 @@ public:
         {
             // element geometry
             const Geometry& geometry = element.geometry();
-
-            Dune::GeometryType geomType = geometry.type();
-
-            const Dune::FieldVector<Scalar,dim>& local = ReferenceElements::general(geomType).position(0, 0);
+            const Dune::FieldVector<Scalar,dim>& local = referenceElement(geometry).position(0, 0);
             Dune::FieldVector<Scalar,dim> globalPos = geometry.global(local);
 
             Scalar volume = geometry.volume();
@@ -709,8 +699,6 @@ public:
         using Element = typename Grid::template Codim<0>::Entity;
         using Geometry = typename Element::Geometry;
         using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
-        using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
-        using ReferenceFaces = Dune::ReferenceElements<Scalar, dim-1>;
 
         ElementMapper elementMapper(gridView, Dune::mcmgElementLayout());
 
@@ -728,14 +716,10 @@ public:
         {
             // element geometry
             const Geometry& geometry = element.geometry();
-
-            Dune::GeometryType geomType = geometry.type();
-
-            const Dune::FieldVector<Scalar,dim>& local = ReferenceElements::general(geomType).position(0, 0);
+            const Dune::FieldVector<Scalar,dim>& local = referenceElement(geometry).position(0, 0);
             Dune::FieldVector<Scalar,dim> globalPos = geometry.global(local);
 
-            Scalar volume = geometry.integrationElement(local)
-                    *ReferenceElements::general(geomType).volume();
+            Scalar volume = geometry.integrationElement(local) * referenceElement(geometry).volume();
 
             int eIdx = elementMapper.index(element);
 
@@ -763,14 +747,11 @@ public:
             Dune::FieldVector<Scalar,dim> exactGradient;
             for (const auto& intersection : intersections(gridView, element))
             {
-                // get geometry type of face
-                Dune::GeometryType gtf = intersection.geometryInInside().type();
-
                 // local number of facet
                 i++;
 
                 // center in face's reference element
-                const Dune::FieldVector<Scalar,dim-1>& faceLocalNm1 = ReferenceFaces::general(gtf).position(0,0);
+                const Dune::FieldVector<Scalar,dim-1>& faceLocalNm1 = referenceElement(intersection.geometryInInside()).position(0,0);
 
                 // center of face in global coordinates
                 Dune::FieldVector<Scalar,dim> faceGlobal = intersection.geometry().global(faceLocalNm1);
-- 
GitLab