add mpfa framework

TODO: implement ficks and fouriers law, as well as the interaction volume data handles
accordingly for problems involving diffusion/heat conduction

dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseeds.hh → dumux/discretization/cellcentered/mpfa/connectivitymap.hh

@@ -18,35 +18,36 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief Base class for the global interaction volume seeds of mpfa methods.
+ * \brief Stores the face indices corresponding to the neighbors of an element
+ *        that contribute to the derivative calculation
  */
-#ifndef DUMUX_DISCRETIZATION_MPFA_GLOBALINTERACTIONVOLUMESEEDS_HH
-#define DUMUX_DISCRETIZATION_MPFA_GLOBALINTERACTIONVOLUMESEEDS_HH
+#ifndef DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
+#define DUMUX_CC_MPFA_CONNECTIVITY_MAP_HH
 
-#include <dumux/implicit/cellcentered/mpfa/properties.hh>
-#include <dumux/discretization/cellcentered/mpfa/globalinteractionvolumeseedsbase.hh>
-#include "methods.hh"
+#include <dumux/discretization/cellcentered/connectivitymap.hh>
+#include <dumux/discretization/cellcentered/mpfa/generalconnectivitymap.hh>
 
 namespace Dumux
 {
-//! forward declaration of the actual method-specific implementation
-//! By default we simply inherit from the base class
-//! Actual implementations for other methods have to be provided below
+//! Forward declaration of method specific implementation of the assembly map
 template<class TypeTag, MpfaMethods method>
-class CCMpfaGlobalInteractionVolumeSeedsImplementation;
+class CCMpfaConnectivityMapImplementation;
 
-/*!
- * \ingroup Mpfa
- * \brief Base class for the creation and storage of the interaction volume seeds for mpfa methods.
- */
+// //! The Assembly map for models using mpfa methods
 template<class TypeTag>
-using CCMpfaGlobalInteractionVolumeSeeds = CCMpfaGlobalInteractionVolumeSeedsImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
+using CCMpfaConnectivityMap = CCMpfaConnectivityMapImplementation<TypeTag, GET_PROP_VALUE(TypeTag, MpfaMethod)>;
+
+//! The default is the general assembly map for mpfa schemes
+template<class TypeTag, MpfaMethods method>
+class CCMpfaConnectivityMapImplementation : public CCMpfaGeneralConnectivityMap<TypeTag> {};
 
-} // end namespace
+//! The o-method can use the simple assembly map
+template<class TypeTag>
+class CCMpfaConnectivityMapImplementation<TypeTag, MpfaMethods::oMethod> : public CCSimpleConnectivityMap<TypeTag> {};
 
-// the specializations of this class differing from the default have to be included here
-#include <dumux/discretization/cellcentered/mpfa/omethod/globalinteractionvolumeseeds.hh>
-#include <dumux/discretization/cellcentered/mpfa/lmethod/globalinteractionvolumeseeds.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethodfps/globalinteractionvolumeseeds.hh>
+//! The o-method with full pressure support can use the simple assembly map
+template<class TypeTag>
+class CCMpfaConnectivityMapImplementation<TypeTag, MpfaMethods::oMethodFps> : public CCSimpleConnectivityMap<TypeTag> {};
+}
 
 #endif

dumux/discretization/cellcentered/mpfa/darcyslaw.hh

@@ -25,8 +25,7 @@
 #ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 #define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 
-#include <dumux/implicit/properties.hh>
-#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
 
 namespace Dumux
 {
@@ -63,63 +62,57 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
 
     // Always use the dynamic type for vectors (compatibility with the boundary)
-    using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-    using CoefficientVector = typename BoundaryInteractionVolume::Vector;
+    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 bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
-    static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
 
     //! 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 BoundaryInteractionVolume::GlobalIndexSet;
-        using PositionVector = typename BoundaryInteractionVolume::PositionVector;
+        using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
+        using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
 
     public:
         //! update cached objects
         template<class InteractionVolume>
-        void updateAdvection(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
+        void updateAdvection(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf)
         {
             const auto& localFaceData = iv.getLocalFaceData(scvf);
-            // update the quantities that are equal for all phases
-            advectionVolVarsStencil_ = iv.volVarsStencil();
-            advectionVolVarsPositions_ = iv.volVarsPositions();
-            advectionTij_ = iv.getTransmissibilities(localFaceData);
 
-            // The neumann fluxes always have to be set per phase
-            for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                phaseNeumannFluxes_[phaseIdx] = iv.getNeumannFlux(localFaceData, phaseIdx);
+            // update the quantities that are equal for all phases
+            advectionSwitchFluxSign_ = localFaceData.isOutside;
+            advectionTij_ = &iv.getTransmissibilities(scvf, localFaceData, dataHandle);
+            advectionVolVarsStencil_ = &dataHandle.volVarsStencil();
+            advectionDirichletData_ = &dataHandle.dirichletData();
         }
 
-        //! Returns the volume variables indices necessary for flux computation
-        //! This includes all participating boundary volume variables. Since we
-        //! do not allow mixed BC for the mpfa this is the same for all phases.
+        //! Returns the stencil for advective scvf flux computation
         const Stencil& advectionVolVarsStencil() const
-        { return advectionVolVarsStencil_; }
-
-        //! Returns the position on which the volume variables live. This is
-        //! necessary as we need to evaluate gravity also for the boundary volvars
-        const PositionVector& advectionVolVarsPositions() const
-        { return advectionVolVarsPositions_; }
+        { return *advectionVolVarsStencil_; }
 
         //! 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_; }
+        { return *advectionTij_; }
 
-        //! If the useTpfaBoundary property is set to false, the boundary conditions
-        //! are put into the local systems leading to possible contributions on all faces
-        Scalar advectionNeumannFlux(unsigned int phaseIdx) const
-        { return phaseNeumannFluxes_[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_; }
+
+        //! Returns the data on dirichlet boundary conditions affecting
+        //! the flux computation on this face
+        const DirichletDataContainer& advectionDirichletData() const
+        { return *advectionDirichletData_; }
 
     private:
-        // Quantities associated with advection
-        Stencil advectionVolVarsStencil_;
-        PositionVector advectionVolVarsPositions_;
-        CoefficientVector advectionTij_;
-        std::array<Scalar, numPhases> phaseNeumannFluxes_;
+        bool advectionSwitchFluxSign_;
+        const Stencil* advectionVolVarsStencil_;
+        const CoefficientVector* advectionTij_;
+        const DirichletDataContainer* advectionDirichletData_;
     };
 
     //! Class that fills the cache corresponding to mpfa Darcy's Law
@@ -138,10 +131,14 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
                          const FluxVariablesCacheFiller& fluxVarsCacheFiller)
         {
             // get interaction volume from the flux vars cache filler & upate the cache
-            if (problem.model().fvGridGeometry().isInBoundaryInteractionVolume(scvf))
-                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf);
+            if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
+                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.secondaryInteractionVolume(),
+                                                  fluxVarsCacheFiller.dataHandle(),
+                                                  scvf);
             else
-                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.interactionVolume(), scvf);
+                scvfFluxVarsCache.updateAdvection(fluxVarsCacheFiller.primaryInteractionVolume(),
+                                                  fluxVarsCacheFiller.dataHandle(),
+                                                  scvf);
         }
     };
 
@@ -164,27 +161,19 @@ public:
         static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
 
         const auto& fluxVarsCache = elemFluxVarsCache[scvf];
-        const auto& volVarsStencil = fluxVarsCache.advectionVolVarsStencil();
-        const auto& volVarsPositions = fluxVarsCache.advectionVolVarsPositions();
         const auto& tij = fluxVarsCache.advectionTij();
 
-        const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
-        // For interior Neumann boundaries when using Tpfa on boundaries, return the user-specified flux
-        // We assume phaseIdx = eqIdx here.
-        if (isInteriorBoundary
-            && useTpfaBoundary
-            && fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
-            return scvf.area()*
-                   elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
-                   problem.neumann(element, fvGeometry, elemVolVars, scvf)[phaseIdx];
-
         // Calculate the interface density for gravity evaluation
-        const auto rho = Implementation::interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
+        const auto rho = interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx);
 
-        // calculate Tij*pj
-        Scalar flux(0.0);
-        unsigned int localIdx = 0;
-        for (const auto volVarIdx : volVarsStencil)
+        // Variable for the flux to be computed
+        Scalar scvfFlux(0.0);
+
+        // index counter to get corresponding transmissibilities
+        unsigned int i = 0;
+
+        // add contributions from cell-centered unknowns
+        for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil())
         {
             const auto& volVars = elemVolVars[volVarIdx];
             Scalar h = volVars.pressure(phaseIdx);
@@ -193,32 +182,41 @@ public:
             if (gravity)
             {
                 // gravitational acceleration in the center of the actual element
-                const auto x = volVarsPositions[localIdx];
+                const auto x = fvGeometry.scv(volVarIdx).center();
                 const auto g = problem.gravityAtPos(x);
-
                 h -= rho*(g*x);
             }
 
-            flux += tij[localIdx++]*h;
+            scvfFlux += tij[i++]*h;
         }
 
-        // if no interior boundaries are present, return the flux
-        if (!enableInteriorBoundaries)
-            return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
+        // 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);
+            }
 
-        // Handle interior boundaries
-        flux += Implementation::computeInteriorBoundaryContribution(problem, fvGeometry, elemVolVars, fluxVarsCache, phaseIdx, rho);
+            scvfFlux += tij[i++]*h;
+        }
 
-        // return overall resulting flux
-        return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
+        // return the flux (maybe adjust the sign)
+        return fluxVarsCache.advectionSwitchFluxSign() ? -scvfFlux : scvfFlux;
     }
 
+private:
     static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
                                      const ElementVolumeVariables& elemVolVars,
                                      const SubControlVolumeFace& scvf,
                                      const FluxVariablesCache& fluxVarsCache,
-                                     const unsigned int phaseIdx,
-                                     const bool isInteriorBoundary)
+                                     const unsigned int phaseIdx)
     {
         static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
 
@@ -226,14 +224,6 @@ public:
             return Scalar(0.0);
         else
         {
-            // Treat interior Dirichlet boundaries differently
-            if (isInteriorBoundary)
-            {
-                const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
-                if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
-                    return data.facetVolVars(fvGeometry).density(phaseIdx);
-            }
-
             // use arithmetic mean of the densities around the scvf
             if (!scvf.boundary())
             {
@@ -246,44 +236,6 @@ public:
                 return elemVolVars[scvf.outsideScvIdx()].density(phaseIdx);
         }
     }
-
-    static Scalar computeInteriorBoundaryContribution(const Problem& problem,
-                                                      const FVElementGeometry& fvGeometry,
-                                                      const ElementVolumeVariables& elemVolVars,
-                                                      const FluxVariablesCache& fluxVarsCache,
-                                                      unsigned int phaseIdx, Scalar rho)
-    {
-        static const bool gravity = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
-
-        // obtain the transmissibilites associated with all pressures
-        const auto& tij = fluxVarsCache.advectionTij();
-
-        // the interior dirichlet boundaries local indices start after
-        // the cell and the domain Dirichlet boundary pressures
-        const auto startIdx = fluxVarsCache.advectionVolVarsStencil().size();
-
-        // add interior Dirichlet boundary contributions
-        Scalar flux = 0.0;
-        for (auto&& data : fluxVarsCache.interiorBoundaryData())
-        {
-            if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
-            {
-                Scalar h = data.facetVolVars(fvGeometry).pressure(phaseIdx);
-
-                if (gravity)
-                {
-                    const auto x = fvGeometry.scvf(data.scvfIndex()).ipGlobal();
-                    const auto g = problem.gravityAtPos(x);
-
-                    h -= rho*(g*x);
-                }
-
-                flux += tij[startIdx + data.localIndexInInteractionVolume()]*h;
-            }
-        }
-
-        return flux;
-    }
 };
 
 } // end namespace

dumux/discretization/cellcentered/mpfa/dualgridindexset.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
+ * \brief Base 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/implicit/cellcentered/mpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Mpfa
+ * \brief Nodal index set for the dual grid of mpfa schemes.
+ */
+template<class TypeTag>
+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;
+
+    static const int dim = GridView::dimension;
+
+public:
+    //! constructor
+    DualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
+
+    // Inserts data for a given scvf. This should only called once per scvf!
+    void insert(const SubControlVolumeFace& scvf)
+    {
+        insert(scvf.boundary(),
+               scvf.index(),
+               scvf.insideScvIdx(),
+               scvf.outsideScvIndices());
+    }
+
+    // Inserts data for a given scvf. This should only called once per scvf!
+    void insert(const bool boundary,
+                const IndexType scvfIdx,
+                const IndexType insideScvIdx,
+                const std::vector<IndexType>& outsideScvIndices)
+    {
+        // check if it this really hasn't been called yet
+        assert( std::find_if( 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;
+        scvIndices.reserve(outsideScvIndices.size()+1);
+        scvIndices.push_back(insideScvIdx);
+        scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
+
+        // if scvf is on boundary, increase counter
+        if (boundary)
+            numBoundaryScvfs_++;
+
+        scvfIndices_.push_back(scvfIdx);
+        scvfIsOnBoundary_.push_back(boundary);
+        scvfNeighborScvIndicesGlobal_.emplace_back( std::move(scvIndices) );
+
+        // 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())
+        {
+            const auto localScvIdx = std::distance(scvIndices_.begin(), it);
+            localScvfIndicesInScv_[localScvIdx].push_back(curScvfLocalIdx);
+        }
+        else
+        {
+            LocalIndexContainer localScvfIndices;
+            localScvfIndices.reserve(dim);
+            localScvfIndices.push_back(curScvfLocalIdx);
+            localScvfIndicesInScv_.emplace_back( std::move(localScvfIndices) );
+            scvIndices_.push_back(insideScvIdx);
+        }
+    }
+
+    // This should be called AFTER (!) all scvf data has been inserted
+    void makeLocal()
+    {
+        // make sure the created index set makes sense so far
+        assert(checkGlobalIndexSet_());
+
+        // compute local neighboring scv indices for the scvfs
+        scvfNeighborScvIndices_.resize(numScvfs());
+        for (unsigned int i = 0; i < numScvfs(); ++i)
+        {
+            const auto& neighborsGlobal = scvfNeighborScvIndicesGlobal_[i];
+            const auto numNeighbors = scvfIsOnBoundary_[i] ? 1 : neighborsGlobal.size();
+
+            scvfNeighborScvIndices_[i].resize(numNeighbors);
+            for (unsigned int j = 0; j < numNeighbors; ++j)
+                scvfNeighborScvIndices_[i][j] = findLocalScvIdx_(neighborsGlobal[j]);
+        }
+
+        // delete global neighboring scv data (not used anymore)
+        scvfNeighborScvIndicesGlobal_.clear();
+    }
+
+    //! returns the number of scvs
+    std::size_t numScvs() const
+    { return scvIndices_.size(); }
+
+    //! returns the number of scvfs
+    std::size_t numScvfs() const
+    { return scvfIndices_.size(); }
+
+    //! returns the number of boundary scvfs
+    std::size_t numBoundaryScvfs() const
+    { return numBoundaryScvfs_; }
+
+    //! returns the global scv indices connected to this dual grid node
+    const GlobalIndexContainer& globalScvIndices() const
+    { return scvIndices_; }
+
+    //! returns the global scvf indices connected to this dual grid node
+    const GlobalIndexContainer& globalScvfIndices() const
+    { return scvfIndices_; }
+
+    //! returns a global scv idx for a given local scv index
+    IndexType scvIdxGlobal(LocalIndexType scvIdxLocal) const
+    { return scvIndices_[scvIdxLocal]; }
+
+    //! returns the local indices of the scvfs embedded in a local scv
+    const LocalIndexContainer& localScvfIndicesInScv(LocalIndexType scvIdxLocal) const
+    { return localScvfIndicesInScv_[scvIdxLocal]; }
+
+    //! returns a global scvf idx for a given local scvf index
+    IndexType scvfIdxGlobal(LocalIndexType scvfIdxLocal) const
+    { return scvfIndices_[scvfIdxLocal]; }
+
+    //! returns whether or not an scvf touches the boundary
+    bool scvfIsOnBoundary(LocalIndexType scvfIdxLocal) const
+    { return scvfIsOnBoundary_[scvfIdxLocal]; }
+
+    //! returns the local indices of the neighboring scvs of an scvf
+    const LocalIndexContainer& localNeighboringScvIndices(LocalIndexType scvfIdxLocal) const
+    { return scvfNeighborScvIndices_[scvfIdxLocal]; }
+
+private:
+    //! returns the node-local scv index to a given global scv index
+    unsigned int findLocalScvIdx_(IndexType globalScvIdx) const
+    {
+        auto it = std::find( scvIndices_.begin(), scvIndices_.end(), globalScvIdx );
+        assert(it != scvIndices_.end() && "Global scv index not found in local container!");
+        return std::distance(scvIndices_.begin(), it);
+    }
+
+    //! checks whether or not the inserted index set makes sense
+    bool checkGlobalIndexSet_() const
+    {
+        //! do all scvs have dim embedded scvfs?
+        for (const auto& scvfIndices : localScvfIndicesInScv_)
+            if (scvfIndices.size() != dim)
+                DUNE_THROW(Dune::InvalidStateException, "Number of scv faces found for an scv != dimension");
+
+        //! are all scvfs unique?
+        for (const auto& scvfIdx : scvfIndices_)
+            for (const auto& scvfIdx2 : scvfIndices_)
+                if (scvfIdx == scvfIdx2)
+                    DUNE_THROW(Dune::InvalidStateException, "Inserted scvfs seem to not be unique");
+
+        return true;
+    }
+
+    //! the indices of the scvs around a dual grid node
+    GlobalIndexContainer scvIndices_;
+    //! maps to each scv a list of scvf indices embedded in it
+    std::vector<LocalIndexContainer> localScvfIndicesInScv_;
+
+    //! the indices of the scvfs around a dual grid node
+    GlobalIndexContainer scvfIndices_;
+    //! maps to each scvf a boolean to indicate if it is on the boundary
+    std::vector<bool> scvfIsOnBoundary_;
+    //! maps to each scvf a list of neighbouring local scv indices
+    //! ordering: 0 - inside scv idx; 1..n - outside scv indices
+    std::vector<LocalIndexContainer> scvfNeighborScvIndices_;
+    //! maps to each scvf a list of neighbouring global scv indices
+    //! This container is destroyed when makeLocal() is called
+    std::vector<GlobalIndexContainer> scvfNeighborScvIndicesGlobal_;
+    //! stores how many boundary scvfs are embedded in this dual grid node
+    std::size_t numBoundaryScvfs_;
+};
+
+/*!
+ * \ingroup Mpfa
+ * \brief Class for the index sets of the dual grid in mpfa schemes.
+ */
+template<class TypeTag>
+class CCMpfaDualGridIndexSet : public std::vector<DualGridNodalIndexSet<TypeTag>>
+{
+    using ParentType = std::vector<DualGridNodalIndexSet<TypeTag>>;
+
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+    static const int dim = GridView::dimension;
+
+public:
+    using NodalIndexSet = DualGridNodalIndexSet<TypeTag>;
+
+    CCMpfaDualGridIndexSet(const GridView& gridView)
+    {
+        (*this).resize(gridView.size(dim));
+    }
+
+    //! after the global data has been inserted, this
+    //! method should be called to set up the local index sets
+    void makeLocalIndexSets()
+    {
+        for (auto& node : (*this))
+            node.makeLocal();
+    }
+
+    //! access with an scvf
+    const NodalIndexSet& operator[] (const SubControlVolumeFace& scvf) const
+    { return (*this)[scvf.vertexIndex()]; }
+    NodalIndexSet& operator[] (const SubControlVolumeFace& scvf)
+    { return (*this)[scvf.vertexIndex()]; }
+
+    //! access with an index
+    using ParentType::operator[];
+};
+} // end namespace
+