[mpfa] restructure framework

Correct handling of gravitational flux contributions requires the
mpfa framework to be completely re-designed.

dumux/discretization/cellcentered/mpfa/CMakeLists.txt

 add_subdirectory("omethod")
 
 install(FILES
+computetransmissibility
 connectivitymap.hh
 darcyslaw.hh
 dualgridindexset.hh
@@ -18,6 +19,8 @@ helper.hh
 interactionvolumebase.hh
 interactionvolumedatahandle.hh
 interactionvolume.hh
+localassembler.hh
+localfacedata.hh
 methods.hh
 properties.hh
 subcontrolvolumeface.hh

dumux/discretization/cellcentered/mpfa/computetransmissibility.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCMpfaDiscretization
+ * \brief This file contains free functions to evaluate the transmissibilities
+ *        associated with flux evaluations across sub-control volume faces
+ *        in the context of cell-centered Mpfa schemes.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_COMPUTE_TRANSMISSIBILITY_HH
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief Free function to evaluate the Mpfa transmissibility associated
+ *        with the flux (in the form of flux = T*gradU) across a
+ *        sub-control volume face stemming from a given sub-control
+ *        volume with corresponding tensor T.
+ *
+ * \param scv The iv-local sub-control volume
+ * \param scvf The grid sub-control volume face
+ * \param T The tensor living in the scv
+ * \param extrusionFactor The extrusion factor of the scv
+ */
+template<class Scv, class Scvf, class Tensor>
+Dune::FieldVector< typename Tensor::field_type, Scv::myDimension >
+computeMpfaTransmissibility(const Scv& scv,
+                            const Scvf& scvf,
+                            const Tensor& T,
+                            typename Scv::ctype extrusionFactor)
+{
+    Dune::FieldVector< typename Tensor::field_type, Scv::myDimension > wijk;
+    for (unsigned int dir = 0; dir < Scv::myDimension; ++dir)
+        wijk[dir] = vtmv(scvf.unitOuterNormal(), T, scv.nu(dir));
+
+    wijk *= scvf.area()*extrusionFactor;
+    wijk /= scv.detX();
+
+    return wijk;
+}
+
+/*!
+ * \ingroup CCMpfaDiscretization
+ * \brief Free function to evaluate the Mpfa transmissibility associated
+ *        with the flux (in the form of flux = T*gradU) across a
+ *        sub-control volume face stemming from a given sub-control
+ *        volume with corresponding tensor T, where T is a scalar.
+ *
+ * \param scv The iv-local sub-control volume
+ * \param scvf The grid sub-control volume face
+ * \param t the scalar quantity living in the scv
+ * \param extrusionFactor The extrusion factor of the scv
+ */
+template< class Scv,
+          class Scvf,
+          class Tensor,
+          std::enable_if_t< Dune::IsNumber<Tensor>::value, int > = 1 >
+Dune::FieldVector<Tensor, Scv::myDimension>
+computeMpfaTransmissibility(const Scv& scv,
+                            const Scvf& scvf,
+                            Tensor t,
+                            typename Scv::ctype extrusionFactor)
+{
+    Dune::FieldVector<Tensor, Scv::myDimension> wijk;
+    for (unsigned int dir = 0; dir < Scv::myDimension; ++dir)
+        wijk[dir] = vtmv(scvf.unitOuterNormal(), t, scv.nu(dir));
+
+    wijk *= scvf.area()*extrusionFactor;
+    wijk /= scv.detX();
+
+    return wijk;
+}
+
+} // end namespace Dumux
+
+#endif

dumux/discretization/cellcentered/mpfa/connectivitymap.hh

@@ -18,8 +18,11 @@
  *****************************************************************************/
 /*!
  * \file
+ * \ingroup CCMpfaDiscretization
  * \brief Stores the face indices corresponding to the neighbors of an element
- *        that contribute to the derivative calculation
+ *        that contribute to the derivative calculation. Depending on if an
+ *        mpfa scheme leads to a symmetric/unsymmetric sparsity pattern, the
+ *        adequate implementation of the connectiviy map is chosen.
  */
 #ifndef DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
 #define DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
@@ -35,15 +38,15 @@ namespace Dumux
 template<class TypeTag, MpfaMethods method>
 class CCMpfaConnectivityMapImplementation;
 
-// //! The Assembly map for models using mpfa methods
+//! The Assembly map for models using mpfa methods
 template<class TypeTag>
 using CCMpfaConnectivityMap = CCMpfaConnectivityMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
 
-//! The default is the general assembly map for mpfa schemes
+//! The default is the general assembly map for mpfa schemes (non-symmetric schemes)
 template<class TypeTag, MpfaMethods method>
 class CCMpfaConnectivityMapImplementation : public CCMpfaGeneralConnectivityMap<TypeTag> {};
 
-//! The o-method can use the simple assembly map
+//! The o-method can use the simple (symmetric) assembly map
 template<class TypeTag>
 class CCMpfaConnectivityMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleConnectivityMap<TypeTag> {};
 }

dumux/discretization/cellcentered/mpfa/darcyslaw.hh

@@ -18,27 +18,29 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief This file contains the class which is required to calculate
- *        volume and mass fluxes of fluid phases over a face of a finite volume by means
- *        of the Darcy approximation. This specializations is for cell-centered schemes
- *        using multi-point flux approximation.
+ * \ingroup CCMpfaDiscretization
+ * \brief Darcy's Law for cell-centered finite volume schemes
+ *        with multi-point flux approximation.
  */
 #ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 #define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 
+#include <dune/common/densevector.hh>
+#include <dune/common/densematrix.hh>
+
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/discretization/methods.hh>
 
 namespace Dumux
 {
-// forward declarations
+//! forward declaration of the method-specific implementation
 template<class TypeTag, DiscretizationMethods discMethod>
 class DarcysLawImplementation;
 
 /*!
- * \ingroup Mpfa
- * \brief Specialization of Darcy's Law for the CCMpfa method.
+ * \ingroup CCMpfaDiscretization
+ * \brief Darcy's law for cell-centered finite volume schemes with multi-point flux approximation.
  */
 template <class TypeTag>
 class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
@@ -47,21 +49,13 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
+
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
 
-    // Always use the dynamic type for vectors (compatibility with the boundary)
-    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
-    using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector;
-    using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
-
-    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
-    static constexpr int dim = GridView::dimension;
-    static constexpr int dimWorld = GridView::dimensionworld;
-
     //! Class that fills the cache corresponding to mpfa Darcy's Law
     class MpfaDarcysLawCacheFiller
     {
@@ -77,14 +71,16 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
                          const SubControlVolumeFace& scvf,
                          const FluxVariablesCacheFiller& fluxVarsCacheFiller)
         {
-            // get interaction volume from the flux vars cache filler & upate the cache
+            // get interaction volume related data from the filler class & upate the cache
             if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
                 scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
-                                                  fluxVarsCacheFiller.dataHandle(),
+                                                  fluxVarsCacheFiller.secondaryIvLocalFaceData(),
+                                                  fluxVarsCacheFiller.secondaryIvDataHandle(),
                                                   scvf);
             else
                 scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
-                                                  fluxVarsCacheFiller.dataHandle(),
+                                                  fluxVarsCacheFiller.primaryIvLocalFaceData(),
+                                                  fluxVarsCacheFiller.primaryIvDataHandle(),
                                                   scvf);
         }
     };
@@ -92,54 +88,124 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
     //! The cache used in conjunction with the mpfa Darcy's Law
     class MpfaDarcysLawCache
     {
-        // We always use the dynamic types here to be compatible on the boundary
-        using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
-        using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
+        static constexpr int dim = GridView::dimension;
+        static constexpr int dimWorld = GridView::dimensionworld;
+        static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+
+        // In the current implementation of the flux variables cache we cannot
+        // make a disctinction between dynamic (mpfa-o) and static (mpfa-l)
+        // matrix and vector types, as currently the cache class can only be templated
+        // by a type tag (and there can only be one). We use a dynamic vector here to
+        // make sure it works in case one of the two used interaction volume types uses
+        // dynamic types performance is thus lowered for schemes using static types.
+        // TODO: this has to be checked thoroughly as soon as a scheme using static types
+        //       is implemented. One idea to overcome the performance drop could be only
+        //       storing the iv-local index here and obtain tij always from the datahandle
+        //       of the fluxVarsCacheContainer
+        using Vector = Dune::DynamicVector< Scalar >;
+        using Matrix = Dune::DynamicMatrix< Scalar >;
+
+        using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+        using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
+        using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
+
+        static_assert( std::is_convertible<PrimaryIvVector*, Vector*>::value,
+                       "The vector type used in primary interaction volumes is not convertible to Dune::DynamicVector!" );
+        static_assert( std::is_convertible<PrimaryIvMatrix*, Matrix*>::value,
+                       "The matrix type used in primary interaction volumes is not convertible to Dune::DynamicMatrix!" );
+
+        using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
+        using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
+        using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
+
+        static_assert( std::is_convertible<SecondaryIvVector*, Vector*>::value,
+                       "The vector type used in secondary interaction volumes is not convertible to Dune::DynamicVector!" );
+        static_assert( std::is_convertible<SecondaryIvMatrix*, Matrix*>::value,
+                       "The matrix type used in secondary interaction volumes is not convertible to Dune::DynamicMatrix!" );
 
     public:
         // export the filler type
         using Filler = MpfaDarcysLawCacheFiller;
 
-        //! update cached objects
-        template<class InteractionVolume>
-        void updateAdvection(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf)
+        /*!
+         * \brief Update cached objects (transmissibilities and gravity)
+         *
+         * \tparam InteractionVolume The (mpfa scheme-specific) interaction volume
+         * \tparam LocalFaceData The object used to store iv-local info on an scvf
+         * \tparam DataHandle The object used to store transmissibility matrices etc.
+         *
+         * \param iv The interaction volume this scvf is embedded in
+         * \param localFaceData iv-local info on this scvf
+         * \param dataHandle Transmissibility matrix & gravity data of this iv
+         * \param scvf The sub-control volume face
+         */
+        template< class InteractionVolume, class LocalFaceData, class DataHandle >
+        void updateAdvection(const InteractionVolume& iv,
+                             const LocalFaceData& localFaceData,
+                             const DataHandle& dataHandle,
+                             const SubControlVolumeFace &scvf)
         {
-            const auto& localFaceData = iv.getLocalFaceData(scvf);
+            switchFluxSign_ = localFaceData.isOutside();
+
+            for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
+                pj_[pIdx] = &dataHandle.pressures(pIdx);
 
-            // update the quantities that are equal for all phases
-            advectionSwitchFluxSign_ = localFaceData.isOutside();
-            advectionVolVarsStencil_ = &dataHandle.volVarsStencil();
-            advectionDirichletData_ = &dataHandle.dirichletData();
+            static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+            const auto ivLocalIdx = localFaceData.ivLocalScvfIndex();
 
-            // the transmissibilities on surface grids have to be obtained from the outside
+            // standard grids
             if (dim == dimWorld)
-                advectionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
+            {
+                Tij_ = &dataHandle.advectionT()[ivLocalIdx];
+
+                if (enableGravity)
+                    for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                        g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
+            }
+
+            // surface grids
             else
-                advectionTij_ = localFaceData.isOutside() ?
-                                &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] :
-                                &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
+            {
+                if (!localFaceData.isOutside())
+                {
+                    Tij_ = &dataHandle.advectionT()[ivLocalIdx];
+
+                    if (enableGravity)
+                        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                            g_[phaseIdx] = dataHandle.gravity(phaseIdx)[ivLocalIdx];
+                }
+                else
+                {
+                    const auto idxInOutsideFaces = localFaceData.scvfLocalOutsideScvfIndex();
+                    Tij_ = &dataHandle.advectionTout()[ivLocalIdx][idxInOutsideFaces];
+
+                    if (enableGravity)
+                        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                            g_[phaseIdx] = dataHandle.gravityOutside(phaseIdx)[ivLocalIdx][idxInOutsideFaces];
+                }
+            }
         }
 
-        //! Returns the stencil for advective scvf flux computation
-        const Stencil& advectionVolVarsStencil() const { return *advectionVolVarsStencil_; }
+        //! Coefficients for the cell (& Dirichlet) unknowns in flux expressions
+        const Vector& advectionTij() const { return *Tij_; }
 
-        //! Returns the transmissibilities associated with the volume variables
-        //! All phases flow through the same rock, thus, tij are equal for all phases
-        const CoefficientVector& advectionTij() const { return *advectionTij_; }
+        //! The cell (& Dirichlet) pressures within this interaction volume
+        const Vector& pressures(unsigned int phaseIdx) const { return *pj_[phaseIdx]; }
 
-        //! On faces that are "outside" w.r.t. a face in the interaction volume,
-        //! we have to take the negative value of the fluxes, i.e. multiply by -1.0
-        bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; }
+        //! The gravitational acceleration for a phase on this scvf
+        Scalar gravity(unsigned int phaseIdx) const { return g_[phaseIdx]; }
 
-        //! Returns the data on dirichlet boundary conditions affecting
-        //! the flux computation on this face
-        const DirichletDataContainer& advectionDirichletData() const { return *advectionDirichletData_; }
+        //! In the interaction volume-local system of eq we have one unknown per face.
+        //! On scvfs on this face, but in "outside" (neighbor) elements of it, we have
+        //! to take the negative value of the fluxes due to the flipped normal vector.
+        //! This function returns whether or not this scvf is an "outside" face in the iv.
+        bool advectionSwitchFluxSign() const { return switchFluxSign_; }
 
     private:
-        bool advectionSwitchFluxSign_;
-        const Stencil* advectionVolVarsStencil_;
-        const CoefficientVector* advectionTij_;
-        const DirichletDataContainer* advectionDirichletData_;
+        bool switchFluxSign_;
+        const Vector* Tij_;                       //!< The transmissibilities such that f = Tij*pj
+        std::array< Scalar, numPhases > g_;       //!< Gravitational flux contribution on this face
+        std::array<const Vector*, numPhases> pj_; //!< The interaction-volume wide phase pressures pj
     };
 
 public:
@@ -149,6 +215,7 @@ public:
     // export the type for the corresponding cache
     using Cache = MpfaDarcysLawCache;
 
+    //! Compute the advective flux across an scvf
     static Scalar flux(const Problem& problem,
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
@@ -157,78 +224,23 @@ public:
                        const unsigned int phaseIdx,
                        const ElementFluxVariablesCache& elemFluxVarsCache)
     {
-        static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
-
-        // Calculate the interface density for gravity evaluation
-        const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
-
-        // prepare computations
-        unsigned int i = 0;
-        Scalar scvfFlux(0.0);
         const auto& fluxVarsCache = elemFluxVarsCache[scvf];
         const auto& tij = fluxVarsCache.advectionTij();
+        const auto& pj = fluxVarsCache.pressures(phaseIdx);
 
-        // add contributions from cell-centered unknowns
-        for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil())
-        {
-            const auto& volVars = elemVolVars[volVarIdx];
-            Scalar h = volVars.pressure(phaseIdx);
+        //! compute t_ij*p_j
+        Scalar scvfFlux = tij*pj;
 
-            // if gravity is enabled, add gravitational acceleration
-            if (gravity)
-            {
-                // gravitational acceleration in the center of the actual element
-                const auto x = fvGeometry.scv(volVarIdx).center();
-                const auto g = problem.gravityAtPos(x);
-                h -= rho*(g*x);
-            }
+        //! maybe add gravitational acceleration
+        static const bool enableGravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
+        if (enableGravity)
+            scvfFlux += fluxVarsCache.gravity(phaseIdx);
 
-            scvfFlux += tij[i++]*h;
-        }
+        //! switch the sign if necessary
+        if (fluxVarsCache.advectionSwitchFluxSign())
+            scvfFlux *= -1.0;
 
-        // add contributions from possible dirichlet boundary conditions
-        for (const auto& d : fluxVarsCache.advectionDirichletData())
-        {
-            const auto& volVars = elemVolVars[d.volVarIndex()];
-            Scalar h = volVars.pressure(phaseIdx);
-
-            // maybe add gravitational acceleration
-            if (gravity)
-            {
-                const auto x = d.ipGlobal();
-                const auto g = problem.gravityAtPos(x);
-                h -= rho*(g*x);
-            }
-
-            scvfFlux += tij[i++]*h;
-        }
-
-        // return the flux (maybe adjust the sign)
-        return fluxVarsCache.advectionSwitchFluxSign() ? -scvfFlux : scvfFlux;
-    }
-
-private:
-    static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
-                                     const SubControlVolumeFace& scvf,
-                                     const unsigned int phaseIdx)
-    {
-        static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
-
-        if (!gravity)
-            return Scalar(0.0);
-        else
-        {
-            // use arithmetic mean of the densities around the scvf
-            if (!scvf.boundary())
-            {
-                Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
-                for (const auto outsideIdx : scvf.outsideScvIndices())
-                    rho += elemVolVars[outsideIdx].density(phaseIdx);
-                return rho/(scvf.outsideScvIndices().size()+1);
-            }
-            else
-                return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
-        }
+        return scvfFlux;
     }
 };
 

dumux/discretization/cellcentered/mpfa/dualgridindexset.hh

@@ -18,39 +18,40 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Base class for the index sets of the dual grid in mpfa schemes.
+ * \ingroup CCMpfaDiscretization
+ * \brief Class for the index sets of the dual grid in mpfa schemes.
  */
 #ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH
 #define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEXSET_BASE_HH
 
-#include <dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh>
-
 namespace Dumux
 {
 
 /*!
- * \ingroup Mpfa
- * \brief Nodal index set for the dual grid of mpfa schemes.
+ * \ingroup CCMpfaDiscretization
+ * \brief Nodal index set for mpfa schemes, constructed
+ *        around grid vertices.
+ *
+ * \tparam GI The type used for indices on the grid
+ * \tparam LI The type used for indexing in interaction volumes
+ * \tparam dim The dimension of the grid
  */
-template<class TypeTag>
+template< class GI, class LI, int dim >
 class DualGridNodalIndexSet
 {
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using IndexType = typename GridView::IndexSet::IndexType;
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
-    using LocalIndexType = typename PrimaryInteractionVolume::Traits::LocalIndexType;
-    using LocalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicLocalIndexContainer;
-    using GlobalIndexContainer = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
+public:
+    using GridIndexType = GI;
+    using LocalIndexType = LI;
 
-    static const int dim = GridView::dimension;
+    // we use dynamic containers to store indices here
+    using GridIndexContainer = std::vector<GridIndexType>;
+    using LocalIndexContainer = std::vector<LocalIndexType>;
 
-public:
-    //! constructor
+    //! Constructor
     DualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
 
-    // Inserts data for a given scvf. This should only called once per scvf!
+    //! Inserts data for a given scvf
+    template<typename SubControlVolumeFace>
     void insert(const SubControlVolumeFace& scvf)
     {
         insert(scvf.boundary(),
@@ -59,20 +60,21 @@ public:
                scvf.outsideScvIndices());
     }
 
-    // Inserts data for a given scvf. This should only called once per scvf!
+    //! Inserts scvf data
     void insert(const bool boundary,
-                const IndexType scvfIdx,
-                const IndexType insideScvIdx,
-                const std::vector<IndexType>& outsideScvIndices)
+                const GridIndexType scvfIdx,
+                const GridIndexType insideScvIdx,
+                const std::vector<GridIndexType>& outsideScvIndices)
     {
-        // this should always be called only once
-        assert( std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() && "scvf has already been inserted!");
+        // this should always be called only once per scvf
+        assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end()
+               && "scvf has already been inserted!");
 
         // the local index of the scvf data about to be inserted
         const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
 
-        // add global scvf data
-        GlobalIndexContainer scvIndices;
+        // add grid scvf data
+        GridIndexContainer scvIndices;
         scvIndices.reserve(outsideScvIndices.size()+1);
         scvIndices.push_back(insideScvIdx);
         scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
@@ -86,7 +88,7 @@ public:
         scvfIsOnBoundary_.push_back(boundary);
         scvfNeighborScvIndices_.emplace_back( std::move(scvIndices) );
 
-        // if entry for the inside scv exists append scvf local index, create otherwise
+        // if entry for the inside scv exists, append scvf local index, create otherwise
         auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
         if (it != scvIndices_.end())
             localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
@@ -100,105 +102,96 @@ public:
         }
     }
 
-    //! returns the number of scvs
+    //! returns the number of scvs around the node
     std::size_t numScvs() const { return scvIndices_.size(); }