diff --git a/dumux/multidomain/facet/box/CMakeLists.txt b/dumux/multidomain/facet/box/CMakeLists.txt
index cf4da19bf2614084d1cfdd4b27b42e93af3a3952..53883f63df589781a5d9ee71a2025c2cc65e63f0 100644
--- a/dumux/multidomain/facet/box/CMakeLists.txt
+++ b/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
diff --git a/dumux/multidomain/facet/box/darcyslaw.hh b/dumux/multidomain/facet/box/darcyslaw.hh
index 672cca361905f0e56c0fa83105e7078f6d304a46..29de16a08888aa1f34d4e6db6349b05ccf4d493b 100644
--- a/dumux/multidomain/facet/box/darcyslaw.hh
+++ b/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)
diff --git a/dumux/multidomain/facet/box/fluxvariablescache.hh b/dumux/multidomain/facet/box/fluxvariablescache.hh
new file mode 100644
index 0000000000000000000000000000000000000000..41593ec3ef31414b4f16e151d3b82d30a7c8b02b
--- /dev/null
+++ b/dumux/multidomain/facet/box/fluxvariablescache.hh
@@ -0,0 +1,113 @@
+// -*- 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
diff --git a/dumux/multidomain/facet/box/properties.hh b/dumux/multidomain/facet/box/properties.hh
index 72303e00bb6da473799360dbf10a9cb13bc2d81d..40e9a7802a7d7a2e19c3059b317d03d55726732d 100644
--- a/dumux/multidomain/facet/box/properties.hh
+++ b/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>