diff --git a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
index 6a2393293c10fc291d13c6351db4dd347f5817fd..c80cc8c692ea5bfd1f6558c58f7f73d7b91a6ffb 100644
--- a/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/darcyslaw.hh
@@ -224,28 +224,9 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = elemVolVars[insideScvIdx];
- // check if we evaluate the permeability in the volume (for discontinuous fields, default)
- // or at the scvf center for analytical permeability fields (e.g. convergence studies)
- using SpatialParams = typename Problem::SpatialParams;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
- auto getPermeability = [&problem](const VolumeVariables& volVars,
- const GlobalPosition& scvfIpGlobal) -> typename SpatialParams::PermeabilityType
- {
- if (SpatialParams::evaluatePermeabilityAtScvfIP())
- {
- // If the permeability at pos method is not overloaded it might have a different return type
- // We do an implicit cast to permeability type in case EvaluatePermeabilityAtScvfIP is not true so that it compiles
- // If it is true make sure that no cast is necessary and the correct return type is specified in the spatial params
- static_assert(!SpatialParams::evaluatePermeabilityAtScvfIP()
- || std::is_same<std::decay_t<typename SpatialParams::PermeabilityType>, std::decay_t<decltype(problem.spatialParams().permeabilityAtPos(scvfIpGlobal))>>::value,
- "permeabilityAtPos doesn't return PermeabilityType stated in the spatial params!");
- return problem.spatialParams().permeabilityAtPos(scvfIpGlobal);
- }
- else
- return volVars.permeability();
- };
-
- const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, getPermeability(insideVolVars, scvf.ipGlobal()),
+
+ const Scalar ti = computeTpfaTransmissibility(scvf, insideScv,
+ getPermeability_(problem, insideVolVars, scvf.ipGlobal()),
insideVolVars.extrusionFactor());
// on the boundary (dirichlet) we only need ti
@@ -260,7 +241,8 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
// refers to the scv of our element, so we use the scv method
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
- const Scalar tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv, getPermeability(outsideVolVars, scvf.ipGlobal()),
+ const Scalar tj = -1.0*computeTpfaTransmissibility(scvf, outsideScv,
+ getPermeability_(problem, outsideVolVars, scvf.ipGlobal()),
outsideVolVars.extrusionFactor());
// harmonic mean (check for division by zero!)
@@ -273,6 +255,21 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
return tij;
}
+
+private:
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return volVars.permeability(); }
+
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
};
/*!
@@ -458,28 +455,9 @@ public:
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = elemVolVars[insideScvIdx];
- // check if we evaluate the permeability in the volume (for discontinuous fields, default)
- // or at the scvf center for analytical permeability fields (e.g. convergence studies)
- using SpatialParams = typename Problem::SpatialParams;
- using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
- auto getPermeability = [&problem](const VolumeVariables& volVars,
- const GlobalPosition& scvfIpGlobal) -> typename SpatialParams::PermeabilityType
- {
- if (SpatialParams::evaluatePermeabilityAtScvfIP())
- {
- // If the permeability at pos method is not overloaded it might have a different return type
- // We do an implicit cast to permeability type in case evaluatePermeabilityAtScvfIP is not true so that it compiles
- // If it is true make sure that no cast is necessary and the correct return type is specified in the spatial params
- static_assert(!SpatialParams::evaluatePermeabilityAtScvfIP()
- || std::is_same<std::decay_t<typename SpatialParams::PermeabilityType>, std::decay_t<decltype(problem.spatialParams().permeabilityAtPos(scvfIpGlobal))>>::value,
- "permeabilityAtPos doesn't return PermeabilityType stated in the spatial params!");
- return problem.spatialParams().permeabilityAtPos(scvfIpGlobal);
- }
- else
- return volVars.permeability();
- };
-
- const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, getPermeability(insideVolVars, scvf.ipGlobal()),
+
+ const Scalar ti = computeTpfaTransmissibility(scvf, insideScv,
+ getPermeability_(problem, insideVolVars, scvf.ipGlobal()),
insideVolVars.extrusionFactor());
// for the boundary (dirichlet) or at branching points we only need ti
@@ -494,7 +472,8 @@ public:
// refers to the scv of our element, so we use the scv method
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
- const Scalar tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability(outsideVolVars, scvf.ipGlobal()),
+ const Scalar tj = computeTpfaTransmissibility(fvGeometry.flipScvf(scvf.index()), outsideScv,
+ getPermeability_(problem, outsideVolVars, scvf.ipGlobal()),
outsideVolVars.extrusionFactor());
// harmonic mean (check for division by zero!)
@@ -507,6 +486,21 @@ public:
return tij;
}
+
+private:
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return volVars.permeability(); }
+
+ template<class Problem, class VolumeVariables,
+ std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0>
+ static decltype(auto) getPermeability_(const Problem& problem,
+ const VolumeVariables& volVars,
+ const GlobalPosition& scvfIpGlobal)
+ { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); }
};
} // end namespace Dumux
diff --git a/dumux/material/spatialparams/fv1p.hh b/dumux/material/spatialparams/fv1p.hh
index 8187c2021b259a74e0f1156fbdc283c2b4e9994e..581970f2fe98dc8979dade635666d4177d79a740 100644
--- a/dumux/material/spatialparams/fv1p.hh
+++ b/dumux/material/spatialparams/fv1p.hh
@@ -28,11 +28,26 @@
#include <dune/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/math.hh>
+#include <dumux/common/typetraits/isvalid.hh>
#include <dune/common/fmatrix.hh>
namespace Dumux {
+#ifndef DOXYGEN
+namespace Detail {
+// helper struct detecting if the user-defined spatial params class has a permeabilityAtPos function
+// for g++ > 5.3, this can be replaced by a lambda
+template<class GlobalPosition>
+struct hasPermeabilityAtPos
+{
+ auto operator()(auto&& a)
+ -> decltype(a.permeabilityAtPos(std::declval<GlobalPosition>()))
+ {};
+};
+} // end namespace Detail
+#endif
+
/*!
* \ingroup SpatialParameters
* \brief The base class for spatial parameters of one-phase problems
@@ -127,24 +142,19 @@ public:
*/
template<class ElementSolution>
decltype(auto) permeability(const Element& element,
- const SubControlVolume& scv,
- const ElementSolution& elemSol) const
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
{
- return asImp_().permeabilityAtPos(scv.center());
- }
+ static_assert(decltype(isValid(Detail::hasPermeabilityAtPos<GlobalPosition>())(this->asImp_()))::value," \n\n"
+ " Your spatial params class has to either implement\n\n"
+ " const PermeabilityType& permeabilityAtPos(const GlobalPosition& globalPos) const\n\n"
+ " or overload this function\n\n"
+ " template<class ElementSolution>\n"
+ " const PermeabilityType& permeability(const Element& element,\n"
+ " const SubControlVolume& scv,\n"
+ " const ElementSolution& elemSol) const\n\n");
- /*!
- * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$
- * \note It is possibly solution dependent.
- *
- * \return permeability
- * \param globalPos The position of the center of the scv
- */
- Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
- {
- DUNE_THROW(Dune::InvalidStateException,
- "The spatial parameters do not provide "
- "a permeability() or permeabilityAtPos() method.");
+ return asImp_().permeabilityAtPos(scv.center());
}
/*!