[facet][box][flux] use custom flux variables cache

This flux variables cache precomputes the shape values at the tpfa
flux support point necessary on interior boundary scvfs with Neumann
BCs specified on them.

dumux/multidomain/facet/box/CMakeLists.txt

@@ -3,6 +3,7 @@ couplingmanager.hh
 couplingmapper.hh
 darcyslaw.hh
 elementboundarytypes.hh
+fluxvariablescache.hh
 fvelementgeometry.hh
 fvgridgeometry.hh
 localresidual.hh

dumux/multidomain/facet/box/darcyslaw.hh

@@ -91,35 +91,15 @@ public:
         // on interior Neumann boundaries, evaluate the flux using the facet permeability
         if (bcTypes.hasOnlyNeumann())
         {
-            // compute point on element whose connection vector to
-            // the scvf integration point is parallel to the normal
-            const auto& elemGeometry = element.geometry();
-            const auto d1 = scvf.ipGlobal() - insideScv.dofPosition();
-            const auto d2 = elemGeometry.center() - insideScv.dofPosition();
-
-            const auto d1Norm = d1.two_norm();
-            const auto d2Norm = d2.two_norm();
-
-            using std::tan;
-            using std::acos;
-            const auto angle = acos( (d1*d2)/d1Norm/d2Norm );
-            const auto dm = tan(angle)*d1Norm;
-
-            auto pos = scvf.unitOuterNormal();
-            pos *= -1.0*dm;
-            pos += scvf.ipGlobal();
-
-            const auto posLocal = elemGeometry.local(pos);
-            std::vector< Dune::FieldVector<Scalar, 1> > insideShapeValues;
-            fvGeometry.fvGridGeometry().feCache().get(elemGeometry.type()).localBasis().evaluateFunction(posLocal, insideShapeValues);
+            const auto& supportPointShapeValues = fluxVarCache.shapeValuesAtTpfaSupportPoint();
 
             Scalar rho = 0.0;
-            Scalar pInside = 0.0;
+            Scalar supportPressure = 0.0;
             for (const auto& scv : scvs(fvGeometry))
             {
                 const auto& volVars = elemVolVars[scv];
                 rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
-                pInside += volVars.pressure(phaseIdx)*insideShapeValues[scv.indexInElement()][0];
+                supportPressure += volVars.pressure(phaseIdx)*supportPointShapeValues[scv.indexInElement()][0];
             }
 
             // compute tpfa flux such that flux and pressure continuity holds
@@ -129,11 +109,12 @@ public:
             using std::sqrt;
             const auto df = dim == dimWorld ? 0.5*facetVolVars.extrusionFactor()
                                             : 0.5*sqrt(facetVolVars.extrusionFactor());
+            const auto dm = (scvf.ipGlobal() - fluxVarCache.tpfaSupportPoint()).two_norm();
 
             const auto tm = 1.0/dm*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), scvf.unitOuterNormal());
             const auto tf = 1.0/df*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal());
 
-            auto flux = tm*tf/(tm+tf)*(pInside - facetVolVars.pressure(phaseIdx))
+            auto flux = tm*tf/(tm+tf)*(supportPressure - facetVolVars.pressure(phaseIdx))
                         *scvf.area()*insideVolVars.extrusionFactor();
 
             if (enableGravity)

dumux/multidomain/facet/box/fluxvariablescache.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 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup BoxDiscretization
+ * \brief Flux variables cache class for the box scheme in the context
+ *        of models considering coupling across the element facets.
+ */
+#ifndef DUMUX_DISCRETIZATION_BOX_FACETCOUPLING_FLUXVARIABLES_CACHE_HH
+#define DUMUX_DISCRETIZATION_BOX_FACETCOUPLING_FLUXVARIABLES_CACHE_HH
+
+#include <cassert>
+
+#include <dune/common/fvector.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dumux/discretization/box/fluxvariablescache.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoxDiscretization
+ * \brief Flux variables cache class for the box scheme in the context
+ *        of models considering coupling across the element facets.
+ *        This class does not contain any physics-/process-dependent
+ *        data. It solely stores disretization-/grid-related data.
+ */
+template< class Scalar, class FVGridGeometry >
+class BoxFacetCouplingFluxVariablesCache
+: public BoxFluxVariablesCache<Scalar, FVGridGeometry>
+{
+    using ParentType = BoxFluxVariablesCache<Scalar, FVGridGeometry>;
+
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+
+    using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
+
+public:
+
+    //! update the cache for an scvf
+    template< class Problem, class ElementVolumeVariables >
+    void update(const Problem& problem,
+                const Element& element,
+                const FVElementGeometry& fvGeometry,
+                const ElementVolumeVariables& elemVolVars,
+                const SubControlVolumeFace& scvf)
+    {
+        ParentType::update(problem, element, fvGeometry, elemVolVars, scvf);
+
+        // on interior boundaries with Neumann BCs, prepare the shape values at a point
+        // inside the element whose orthogonal projection is the integration point on scvf
+        if (scvf.interiorBoundary() && problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann())
+        {
+            isInteriorNeumannCache_ = true;
+            const auto& geometry = element.geometry();
+            const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+
+            const auto d1 = scvf.ipGlobal() - insideScv.dofPosition();
+            const auto d2 = geometry.center() - insideScv.dofPosition();
+
+            const auto d1Norm = d1.two_norm();
+            const auto d2Norm = d2.two_norm();
+
+            using std::tan;
+            using std::acos;
+            const auto angle = acos( (d1*d2)/d1Norm/d2Norm );
+            const auto dm = tan(angle)*d1Norm;
+
+            ipGlobalInside_ = scvf.unitOuterNormal();
+            ipGlobalInside_ *= -1.0*dm;
+            ipGlobalInside_ += scvf.ipGlobal();
+
+            fvGeometry.feLocalBasis().evaluateFunction(geometry.local(ipGlobalInside_), shapeValuesInside_);
+        }
+    }
+
+    //! returns the integration point inside the element for interior boundaries
+    const GlobalPosition& tpfaSupportPoint() const
+    { assert(isInteriorNeumannCache_); return ipGlobalInside_; }
+
+    //! returns the shape values at ip inside the element for interior boundaries
+    const std::vector<ShapeValue>& shapeValuesAtTpfaSupportPoint() const
+    { assert(isInteriorNeumannCache_); return shapeValuesInside_; }
+
+private:
+    bool isInteriorNeumannCache_{false};
+    GlobalPosition ipGlobalInside_;
+    std::vector<ShapeValue> shapeValuesInside_;
+};
+
+} // end namespace Dumux
+
+#endif // DUMUX_DISCRETIZATION_BOX_FACETCOUPLING_FLUXVARIABLES_CACHE_HH

dumux/multidomain/facet/box/properties.hh

@@ -33,6 +33,7 @@
 #include <dumux/discretization/box.hh>
 
 #include <dumux/multidomain/facet/box/darcyslaw.hh>
+#include <dumux/multidomain/facet/box/fluxvariablescache.hh>
 #include <dumux/multidomain/facet/box/elementboundarytypes.hh>
 #include <dumux/multidomain/facet/box/fvgridgeometry.hh>
 #include <dumux/multidomain/facet/box/localresidual.hh>
@@ -63,6 +64,14 @@ struct AdvectionType<TypeTag, TTag::BoxFacetCouplingModel>
                                             GetPropType<TypeTag, Properties::FVGridGeometry> >;
 };
 
+//! Use the flux variables cache specific for box facet coupling models
+template<class TypeTag>
+struct FluxVariablesCache<TypeTag, TTag::BoxFacetCouplingModel>
+{
+    using type = BoxFacetCouplingFluxVariablesCache< GetPropType<TypeTag, Properties::Scalar>,
+                                                     GetPropType<TypeTag, Properties::FVGridGeometry> >;
+};
+
 //! Per default, use the porous medium flow flux variables with the modified upwind scheme
 template<class TypeTag>
 struct FluxVariables<TypeTag, TTag::BoxFacetCouplingModel>