[mpfa] implement the fps scheme in 2d

dumux/common/math.hh

@@ -512,6 +512,22 @@ Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
                      const Dune::FieldVector<Scalar, 3> &vec3)
 {   return crossProduct<Scalar>(vec1, vec2)*vec3; }
 
+/*!
+ * \brief Transpose a FieldMatrix
+ *
+ * \param M The matrix to be transposed
+ */
+template <class Scalar, int m, int n>
+Dune::FieldMatrix<Scalar, m, n> getTransposed(const Dune::FieldMatrix<Scalar, m, n>& M)
+{
+    Dune::FieldMatrix<Scalar, m, n> T;
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < m; ++j)
+            T[j][i] = M[i][j];
+
+    return T;
+}
+
 /*!
  * \brief Multiply two dynamic matrices
  *

dumux/discretization/cellcentered/mpfa/helper.hh

@@ -267,5 +267,6 @@ public:
 
 // The implemented helper classes need to be included here
 #include <dumux/discretization/cellcentered/mpfa/omethod/helper.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh>
 
 #endif

dumux/discretization/cellcentered/mpfa/interactionvolume.hh

@@ -43,5 +43,6 @@ using CCMpfaInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, G
 
 // the specializations of this class for the available methods have to be included here
 #include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
+#include <dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh>
 
 #endif

dumux/discretization/cellcentered/mpfa/methods.hh

@@ -27,7 +27,8 @@ namespace Dumux
 {
     enum class MpfaMethods : unsigned int
     {
-        oMethod
+        oMethod,
+        oMethodFps
     };
 } // end namespace
 

dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh

@@ -104,6 +104,9 @@ public:
 template<class TypeTag, class Traits>
 class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase<TypeTag, Traits>
 {
+    // The interaction volume of the mpfa-o fps method has to be friend
+    friend CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps>;
+
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -337,7 +340,7 @@ private:
         for (const auto& scvSeed : seed.scvSeeds())
         {
             auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
-            localScvs_.emplace_back(LocalScvType(element, fvGeometry_(), scvSeed));
+            localScvs_.emplace_back(LocalScvType(problem_(), element, fvGeometry_(), scvSeed));
             localElements_.emplace_back(std::move(element));
         }
 

dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh

@@ -32,6 +32,7 @@ template<class TypeTag>
 class CCMpfaOLocalScv
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
@@ -51,7 +52,8 @@ class CCMpfaOLocalScv
 
 public:
     // constructor has the same signature as the LocalScv entity
-    CCMpfaOLocalScv(const Element& element,
+    CCMpfaOLocalScv(const Problem& problem,
+                    const Element& element,
                     const FVElementGeometry& fvGeometry,
                     const LocalScvSeed& scvSeed)
     : seed_(scvSeed)

dumux/discretization/cellcentered/mpfa/omethodfps/helper.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 Helper class to get the required information on an interaction volume.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFAO_FPS_HELPER_HH
+#define DUMUX_DISCRETIZATION_CC_MPFAO_FPS_HELPER_HH
+
+#include <dumux/discretization/cellcentered/mpfa/omethod/helper.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Mpfa
+ * \brief Helper class to get the required information on an interaction volume.
+ *        Specialization for the Mpfa-O method in two dimensions.
+ */
+template<class TypeTag, int dim>
+class MpfaHelperBase<TypeTag, MpfaMethods::oMethodFps, dim> : public MpfaHelperBase<TypeTag, MpfaMethods::oMethod, dim>
+{};
+
+} // end namespace
+
+#endif

dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.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 classes for interaction volume of mpfa methods.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_INTERACTIONVOLUME_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_INTERACTIONVOLUME_HH
+
+#include <dumux/common/math.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh>
+#include <dumux/discretization/cellcentered/mpfa/interactionvolumeseed.hh>
+#include <dumux/discretization/cellcentered/mpfa/facetypes.hh>
+#include <dumux/discretization/cellcentered/mpfa/methods.hh>
+
+#include "localsubcontrolentities.hh"
+
+namespace Dumux
+{
+//! Specialization of the interaction volume traits class
+template<class TypeTag>
+class CCMpfaOFpsInteractionVolumeTraits : public CCMpfaOInteractionVolumeTraits<TypeTag>
+{
+public:
+    // the fps method uses its own interaction volumes at the boundary
+    using BoundaryInteractionVolume = CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps>;
+
+    // The local sub-control volume type differ from the standard mpfa-o method
+    using LocalScvType = CCMpfaOFpsLocalScv<TypeTag>;
+};
+/*!
+ * \ingroup Mpfa
+ * \brief Base class for the interaction volumes of the mpfa-o method with full pressure support.
+ */
+template<class TypeTag>
+class CCMpfaInteractionVolumeImplementation<TypeTag, MpfaMethods::oMethodFps> : public CCMpfaOInteractionVolume<TypeTag, CCMpfaOFpsInteractionVolumeTraits<TypeTag>>
+{
+    using Traits = CCMpfaOFpsInteractionVolumeTraits<TypeTag>;
+    using ParentType = CCMpfaOInteractionVolume<TypeTag, Traits>;
+
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+    using LocalScvType = typename Traits::LocalScvType;
+    using LocalScvfType = typename Traits::LocalScvfType;
+
+    static const int dim = GridView::dimension;
+    using CoordScalar = typename GridView::ctype;
+    using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
+    using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
+    using ShapeJacobian = typename FeLocalBasis::Traits::JacobianType;
+    using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
+    using JacobianInverseTransposed = typename LocalScvType::Geometry::JacobianInverseTransposed;
+
+    using GlobalIdxType = typename Traits::GlobalIndexType;
+    using LocalIdxType = typename Traits::LocalIndexType;
+    using GlobalPosition = typename Traits::GlobalPosition;
+    using LocalPosition = typename LocalScvType::Geometry::LocalCoordinate;
+    using PosVector = typename Traits::PositionVector;
+    using DynamicVector = typename Traits::Vector;
+    using DynamicMatrix = typename Traits::Matrix;
+    using Tensor = typename Traits::Tensor;
+
+    using IVSeed = typename Traits::Seed;
+
+    enum
+    {
+        jacRows = JacobianInverseTransposed::rows,
+        jacCols = JacobianInverseTransposed::cols
+    };
+
+    // stores the matrices required to calculate Tij
+    struct LocalMatrixContainer
+    {
+        DynamicMatrix AL; //! coefficients of unknown face pressures (LHS)
+        DynamicMatrix AR; //! coefficients of unknown face pressures (RHS)
+        DynamicMatrix BL; //! coefficients of cell/Dirichlet pressures (LHS)
+        DynamicMatrix BR; //! coefficients of cell/Dirichlet pressures (RHS)
+
+        // the matrices for the expression of the fluxes
+        DynamicMatrix AF; //! coefficients for the unknown face pressures
+        DynamicMatrix BF; //! coefficients for the cell/Dirichlet pressures
+    };
+
+public:
+
+    CCMpfaInteractionVolumeImplementation(const IVSeed& seed,
+                                          const Problem& problem,
+                                          const FVElementGeometry& fvGeometry,
+                                          const ElementVolumeVariables& elemVolVars)
+    : ParentType(seed, problem, fvGeometry, elemVolVars),
+      c_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, C)),
+      p_(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Mpfa, P)),
+      divEqIdx_(this->fluxScvfIndexSet_().size())
+    {}
+
+    template<typename GetTensorFunction>
+    void solveLocalSystem(const GetTensorFunction& getTensor)
+    {
+        const std::size_t numFluxFaces = this->fluxScvfIndexSet_().size();
+        const std::size_t numUnknowns = numFluxFaces + 1;
+        const std::size_t numFaces = this->localScvfs_.size();
+        const std::size_t numPotentials = this->volVarsStencil().size();
+
+        // instantiate and resize the local matrices
+        LocalMatrixContainer mc;
+        mc.AL.resize(numUnknowns, numUnknowns, 0.0);
+        mc.AR.resize(numUnknowns, numUnknowns, 0.0);
+        mc.BL.resize(numUnknowns, numPotentials, 0.0);
+        mc.BR.resize(numUnknowns, numPotentials, 0.0);
+        mc.AF.resize(numFaces, numUnknowns, 0.0);
+        mc.BF.resize(numFaces, numPotentials, 0.0);
+
+        // assemble the local eq system and matrices
+        assembleLocalMatrices_(getTensor, mc);
+
+        // solve local system and store transmissibility matrix
+        mc.AL -= mc.AR;
+        mc.BR -= mc.BL;
+        mc.AL.invert();
+        this->T_ = Dumux::multiplyMatrices(mc.AF, Dumux::multiplyMatrices(mc.AL, mc.BR));
+        this->T_ += mc.BF;
+    }
+
+    Scalar getNeumannFlux(const std::pair<LocalIdxType, bool>& localIndexPair) const
+    { return 0.0; }
+
+private:
+
+    template<typename GetTensorFunction>
+    void assembleLocalMatrices_(const GetTensorFunction& getTensor,
+                                LocalMatrixContainer& mc)
+    {
+        // loop over the local faces
+        LocalIdxType localFaceIdx = 0;
+        for (const auto& localScvf : this->localScvfs_)
+        {
+            if (localScvf.boundary())
+                assemblePositiveScv(getTensor, localScvf, localFaceIdx, mc, true);
+            else
+            {
+                assemblePositiveScv(getTensor, localScvf, localFaceIdx, mc);
+                assembleNegativeScv(getTensor, localScvf, localFaceIdx, mc);
+            }
+
+            // go to the next face
+            localFaceIdx++;
+        }
+    }
+
+    template<typename GetTensorFunction>
+    void assemblePositiveScv(const GetTensorFunction& getTensor,
+                             const LocalScvfType& localScvf,
+                             const LocalIdxType localScvfIdx,
+                             LocalMatrixContainer& mc,
+                             const bool boundary = false)
+    {
+        // get diffusion tensor in "positive" sub volume
+        const auto localScvIdx = localScvf.insideLocalScvIndex();
+        const auto& localScv = this->localScv_(localScvIdx);
+        const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
+        const auto& element = this->localElement_(localScvIdx);
+        auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
+
+        // the local finite element basis
+        const auto& localBasis = feCache_.get(localScv.geometry().type()).localBasis();
+
+        // the normal and local integration point
+        // On the ref element, normal vector points in the direction of local coordinate
+        auto normalDir = localScv.getScvfIdxInScv(localScvfIdx);
+        auto ipLocal = localScv.geometry().local(localScvf.ip());
+
+        // find normal coordinate direction and integration point for divergence condition
+        LocalIdxType divEqNormalDir = normalDir == 1 ? 0 : 1;
+        LocalPosition divEqIpLocal(0.0);
+        divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
+        divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
+
+        // does the face has an unknown associated with it?
+        bool isFluxFace = localScvf.faceType() != MpfaFaceTypes::dirichlet && localScvf.faceType() != MpfaFaceTypes::interiorDirichlet;
+
+        // assemble coefficients for the face fluxes
+        addFaceFluxCoefficients_(localScv, localBasis, D, localScvfIdx, ipLocal, normalDir, mc, isFluxFace);
+
+        // assemble matrix entries for the condition of zero divergence
+        addDivEquationCoefficients_(localScv, localBasis, D, divEqIpLocal, divEqNormalDir, mc);
+
+        // on boundary faces, add the fluxes from dirichlet boundary conditions
+        if (boundary && !isFluxFace)
+        {
+            LocalPosition bcIpLocal(0.0);
+            bcIpLocal[normalDir] = 1.0;
+            bcIpLocal[divEqNormalDir] = normalDir == 1 ? 1.0 - (1.0-c_)*p_ : c_ + (1.0-c_)*p_;
+            addDivEquationCoefficients_(localScv, localBasis, D, bcIpLocal, normalDir, mc, boundary);
+        }
+    }
+
+    template<typename GetTensorFunction>
+    void assembleNegativeScv(const GetTensorFunction& getTensor,
+                             const LocalScvfType& localScvf,
+                             const LocalIdxType localScvfIdx,
+                             LocalMatrixContainer& mc)
+    {
+        // get diffusion tensor in "negative" sub volume
+        const auto localScvIdx = localScvf.outsideLocalScvIndex();
+        const auto& localScv = this->localScv_(localScvIdx);
+        const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
+        const auto& element = this->localElement_(localScvIdx);;
+        auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
+
+        // the local finite element bases of the scvs
+        const auto& localBasis = feCache_.get(localScv.geometry().type()).localBasis();
+
+        // the normal and local integration point
+        // On the ref element, normal vector points in the direction of local coordinate
+        auto normalDir = localScv.getScvfIdxInScv(localScvfIdx);
+        auto ipLocal = localScv.geometry().local(localScvf.ip());
+
+        // find normals and integration points in the two scvs for condition of zero divergence
+        LocalIdxType divEqNormalDir = normalDir == 1 ? 0 : 1;
+        LocalPosition divEqIpLocal(0.0);
+        divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
+        divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
+
+        // does the face has an unknown associated with it?
+        bool isFluxFace = localScvf.faceType() != MpfaFaceTypes::dirichlet && localScvf.faceType() != MpfaFaceTypes::interiorDirichlet;
+
+        // assemble coefficients for the face fluxes
+        addFaceFluxCoefficients_(localScv, localBasis, D, localScvfIdx, ipLocal, normalDir, mc, isFluxFace, true);
+
+        // assemble matrix entries for the condition of zero divergence
+        addDivEquationCoefficients_(localScv, localBasis, D, divEqIpLocal, divEqNormalDir, mc);
+    }
+
+    void addFaceFluxCoefficients_(const LocalScvType& localScv,
+                                  const FeLocalBasis& localBasis,
+                                  const Tensor& D,
+                                  const LocalIdxType rowIdx,
+                                  const LocalPosition& ipLocal,
+                                  const LocalIdxType normalDir,
+                                  LocalMatrixContainer& mc,
+                                  const bool isFluxEq,
+                                  const bool isRHS = false)
+    {
+        // In case we're on a flux continuity face, get local index
+        LocalIdxType eqSystemIdx = isFluxEq ? this->findLocalIndex(this->fluxScvfIndexSet_(), rowIdx) : -1;
+
+        // Fluxes stemming from the RHS have to have the opposite sign
+        Scalar factor = isRHS ? 1.0 : -1.0;
+        DynamicMatrix& A = isRHS ? mc.AR : mc.AL;
+        DynamicMatrix& B = isRHS ? mc.BR : mc.BL;
+
+        // evaluate shape functions gradients at the ip
+        std::vector<ShapeJacobian> shapeJacobian;
+        localBasis.evaluateJacobian(ipLocal, shapeJacobian);
+
+        // get the inverse jacobian and its transposed
+        const auto jacInvT = localScv.geometry().jacobianInverseTransposed(ipLocal);
+        const auto jacInv = Dumux::getTransposed(jacInvT);
+
+        // transform the diffusion tensor into local coordinates
+        auto DJinvT = jacInvT.leftmultiplyany(D);
+        auto localD = DJinvT.leftmultiplyany(jacInv);
+        localD *= localScv.geometry().integrationElement(ipLocal);
+
+        // add matrix entries for the pressure in the cell center
+        auto cellPressureIdx = this->findLocalIndex(this->volVarsStencil(), localScv.globalIndex());
+        Scalar bi0 = factor*(localD[normalDir]*shapeJacobian[0][0]);
+        if (!isRHS) mc.BF[rowIdx][cellPressureIdx] += bi0;
+        if (isFluxEq) B[eqSystemIdx][cellPressureIdx] += bi0;
+
+        // Add entries from the local scv faces
+        for (int localDir = 0; localDir < dim; localDir++)
+        {
+            auto localScvfIdx = localScv.localScvfIndex(localDir);
+            auto faceType = this->localScvf_(localScvfIdx).faceType();
+            bool otherHasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
+
+            if (otherHasUnknown)
+            {
+                Scalar aij = factor*(localD[normalDir]*shapeJacobian[localDir+1][0]);
+                auto colIdx = this->findLocalIndex(this->fluxScvfIndexSet_(), localScvfIdx);
+                if (!isRHS) mc.AF[rowIdx][colIdx] += aij;
+                if (isFluxEq) A[eqSystemIdx][colIdx] += aij;
+            }
+            else
+            {
+                Scalar bij = factor*(localD[normalDir]*shapeJacobian[localDir+1][0]);
+                auto colIdx = this->localScvs_.size() + this->findLocalIndex(this->dirichletScvfIndexSet_(), localScvfIdx);
+                if (!isRHS) mc.BF[rowIdx][colIdx] += bij;
+                if (isFluxEq) B[eqSystemIdx][colIdx] += bij;
+            }
+        }
+
+        // add entry from the vertex pressure
+        Scalar ain = factor*(localD[normalDir]*shapeJacobian[3][0]);
+        if (!isRHS) mc.AF[rowIdx][divEqIdx_] += ain;
+        if (isFluxEq) A[eqSystemIdx][divEqIdx_] += ain;
+    }
+
+    void addDivEquationCoefficients_(const LocalScvType& localScv,
+                                     const FeLocalBasis& localBasis,
+                                     const Tensor& D,
+                                     const LocalPosition& ipLocal,
+                                     const LocalIdxType normalDir,
+                                     LocalMatrixContainer& mc,
+                                     const bool isBoundary = false)
+    {
+        // fluxes on the auxiliary volume have to be scaled
+        static const Scalar auxArea = 1.0 - c_;
+        Scalar factor = isBoundary ? -1.0*auxArea : auxArea;
+
+        // evaluate shape functions gradients at the ip
+        std::vector<ShapeJacobian> shapeJacobian;
+        localBasis.evaluateJacobian(ipLocal, shapeJacobian);
+
+        // get the inverse jacobian and its transposed
+        const auto jacInvT = localScv.geometry().jacobianInverseTransposed(ipLocal);
+        const auto jacInv = Dumux::getTransposed(jacInvT);
+
+        // transform the diffusion tensor into local coordinates
+        auto DJinvT = jacInvT.leftmultiplyany(D);
+        auto localD = DJinvT.leftmultiplyany(jacInv);
+        localD *= localScv.geometry().integrationElement(ipLocal);
+
+        // add matrix entries for the pressure in the cell center
+        auto cellPressureIdx = this->findLocalIndex(this->volVarsStencil(), localScv.globalIndex());
+        mc.BL[divEqIdx_][cellPressureIdx] += factor*(localD[normalDir]*shapeJacobian[0][0]);
+
+        // Add entries from the local scv faces
+        for (int localDir = 0; localDir < dim; localDir++)
+        {
+            auto localScvfIdx = localScv.localScvfIndex(localDir);
+            auto faceType = this->localScvf_(localScvfIdx).faceType();
+            bool otherHasUnknown = faceType != MpfaFaceTypes::dirichlet && faceType != MpfaFaceTypes::interiorDirichlet;
+
+            if (otherHasUnknown)
+            {
+                auto colIdx = this->findLocalIndex(this->fluxScvfIndexSet_(), localScvfIdx);
+                mc.AL[divEqIdx_][colIdx] += factor*(localD[normalDir]*shapeJacobian[localDir+1][0]);
+            }
+            else
+            {
+                auto colIdx = this->localScvs_.size() + this->findLocalIndex(this->dirichletScvfIndexSet_(), localScvfIdx);
+                mc.BL[divEqIdx_][colIdx] += factor*(localD[normalDir]*shapeJacobian[localDir+1][0]);
+            }
+        }
+
+        // add entry from the vertex pressure
+        mc.AL[divEqIdx_][divEqIdx_] += factor*(localD[normalDir]*shapeJacobian[3][0]);
+    }
+
+    Tensor makeTensor_(Tensor&& tensor) const
+    { return std::move(tensor); }
+
+    Tensor makeTensor_(Scalar&& t) const
+    { return Tensor(t); }
+
+    const FeCache feCache_;
+
+    const Scalar c_;
+    const Scalar p_;
+    const LocalIdxType divEqIdx_;
+};
+
+} // end namespace
+
+#endif

dumux/discretization/cellcentered/mpfa/omethodfps/localsubcontrolentities.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 sub control entities of the mpfa-o method.
+ */
+#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_LOCALSUBCONTROLENTITIES_HH
+#define DUMUX_DISCRETIZATION_CC_MPFA_O_FPS_LOCALSUBCONTROLENTITIES_HH
+
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
+
+#include "mpfafpsgeometryhelper.hh"
+
+namespace Dumux
+{
+template<class TypeTag>
+class CCMpfaOFpsLocalScv
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Helper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using LocalScvSeed = typename InteractionVolume::Seed::LocalScvSeed;
+
+    using GlobalIndexType = typename GridView::IndexSet::IndexType;
+    using LocalIndexType = typename InteractionVolume::LocalIndexType;
+