From eb4dd0d691123b200d8c7b1b5aef6a95e222508c Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 26 Oct 2020 08:55:51 +0100
Subject: [PATCH] [2p][sequential] Adapt headers

---
 dumux/material/spatialparams/sequentialfv.hh  | 20 +++++---
 .../diffusion/cellcentered/pressure.hh        | 36 +++++++++-----
 .../diffusion/cellcentered/velocity.hh        | 21 ++++++---
 .../sequential/diffusion/mimetic/mimetic.hh   | 14 ++++--
 .../diffusion/mimetic/mimeticadaptive.hh      | 14 ++++--
 .../sequential/diffusion/mimetic/operator.hh  | 14 ++++--
 .../diffusion/mimetic/operatoradaptive.hh     | 14 ++++--
 .../sequential/diffusion/mimetic/pressure.hh  | 16 ++++---
 .../diffusion/mimetic/pressureadaptive.hh     | 16 ++++---
 .../diffusion/mpfa/lmethod/2dpressure.hh      | 31 +++++++-----
 .../mpfa/lmethod/2dpressureadaptive.hh        | 33 +++++++------
 .../mpfa/lmethod/2dpressurevelocity.hh        | 19 +++++---
 .../lmethod/2dpressurevelocityadaptive.hh     | 19 +++++---
 .../diffusion/mpfa/lmethod/2dvelocity.hh      | 20 ++++----
 .../mpfa/lmethod/2dvelocityadaptive.hh        |  1 -
 .../diffusion/mpfa/lmethod/3dpressure.hh      | 28 +++++++----
 .../mpfa/lmethod/3dpressureadaptive.hh        |  1 -
 .../mpfa/lmethod/3dpressurevelocity.hh        | 19 +++++---
 .../lmethod/3dpressurevelocityadaptive.hh     | 19 +++++---
 .../diffusion/mpfa/lmethod/3dvelocity.hh      | 18 +++----
 .../mpfa/lmethod/3dvelocityadaptive.hh        |  1 -
 .../diffusion/mpfa/omethod/2dpressure.hh      | 33 +++++++------
 .../mpfa/omethod/2dpressurevelocity.hh        | 19 +++++---
 .../diffusion/mpfa/omethod/2dvelocity.hh      | 20 ++++----
 .../cellcentered/capillarydiffusion.hh        | 25 +++++++---
 .../cellcentered/evalcflfluxcoats.hh          | 47 ++++++++++++-------
 .../cellcentered/evalcflfluxdefault.hh        | 15 ++++--
 .../transport/cellcentered/gravitypart.hh     | 25 +++++++---
 .../transport/cellcentered/saturation.hh      | 29 ++++++++----
 29 files changed, 368 insertions(+), 219 deletions(-)

diff --git a/dumux/material/spatialparams/sequentialfv.hh b/dumux/material/spatialparams/sequentialfv.hh
index d491465768..4dedf42b25 100644
--- a/dumux/material/spatialparams/sequentialfv.hh
+++ b/dumux/material/spatialparams/sequentialfv.hh
@@ -53,9 +53,6 @@ class SequentialFVSpatialParams: public SequentialFVSpatialParamsOneP<TypeTag>
 
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    /// @cond false
-    using MaterialLawParams = typename GetPropType<TypeTag, Properties::MaterialLaw>::Params;
-    /// @endcond
 
 public:
     SequentialFVSpatialParams(const Problem& problem)
@@ -63,15 +60,26 @@ public:
     {
     }
 
+
+    // make sure we get a deprecation warning (remove this after release 3.3)
+    template<class ElementSolution, class SubControlVolume>
+    [[deprecated("Use the new style material laws. Old material laws and this interface will no longer be supported after release 3.3")]]
+    decltype(auto) materialLawParamsDeprecated(const Element& element,
+                                               const SubControlVolume& scv,
+                                               const ElementSolution& elemSol) const
+    {
+        return this->asImp_().materialLawParams(element);
+    }
+
     /*!
      * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
      *
      * \return the material parameters object
      * \param element The element
      */
-    const MaterialLawParams& materialLawParams(const Element &element) const
+    const auto& materialLawParams(const Element &element) const
     {
-            return asImp_().materialLawParamsAtPos(element.geometry().center());
+        return asImp_().materialLawParamsAtPos(element.geometry().center());
     }
 
     /*!
@@ -80,7 +88,7 @@ public:
      * \return the material parameters object
      * \param globalPos The position of the center of the element
      */
-    const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+    const auto& materialLawParamsAtPos(const GlobalPosition& globalPos) const
     {
         DUNE_THROW(Dune::InvalidStateException,
                    "The spatial parameters do not provide "
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
index b3c867d583..d68266be0d 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -121,7 +123,6 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -920,13 +921,20 @@ const Intersection& intersection, const CellData& cellData, const bool first)
             satW = cellData.saturation(wPhaseIdx);
             satNw = cellData.saturation(nPhaseIdx);
         }
-        Scalar temperature = problem_.temperature(element);
+
+        const Scalar temperature = problem_.temperature(element);
 
         // get Dirichlet pressure boundary condition
-        Scalar pressBound = boundValues[pressureIdx];
+        const Scalar pressBound = boundValues[pressureIdx];
 
         //calculate constitutive relations depending on the kind of saturation used
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressW = 0;
@@ -967,9 +975,9 @@ const Intersection& intersection, const CellData& cellData, const bool first)
             rhoMeanNw = 0.5 * (cellData.density(nPhaseIdx) + densityNwBound);
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaWBound = fluidMatrixInteraction.krw(satW)
                 / viscosityWBound;
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW)
                 / viscosityNwBound;
 
         Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNwBound);
@@ -1095,13 +1103,17 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar temperature = problem_.temperature(element);
+        const Scalar temperature = problem_.temperature(element);
 
         // determine phase saturations from primary saturation variable
+        const Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressW = 0;
@@ -1151,8 +1163,8 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
         }
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[nPhaseIdx];
+        Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         if (compressibility_)
         {
@@ -1168,7 +1180,7 @@ void FVPressure2P<TypeTag>::updateMaterialLaws()
         cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
         cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
 
-        Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
+        const Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
 
         Scalar potW = pressW + gravityDiff * density_[wPhaseIdx];
         Scalar potNw = pressNw + gravityDiff * density_[nPhaseIdx];
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
index 7bd5b3f8a9..337358c6c6 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/velocity.hh
@@ -28,6 +28,8 @@
 #include <dune/grid/common/gridenums.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -64,7 +66,6 @@ class FVVelocity2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -564,8 +565,14 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
             satNw = cellData.saturation(nPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         // determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -582,7 +589,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
         }
 
         // get temperature at current position
-        Scalar temperature = problem_.temperature(element);
+        const Scalar temperature = problem_.temperature(element);
 
         Scalar densityWBound = density_[wPhaseIdx];
         Scalar densityNwBound = density_[nPhaseIdx];
@@ -604,9 +611,9 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
             viscosityNwBound = FluidSystem::viscosity(fluidState, nPhaseIdx) / densityNwBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaWBound = fluidMatrixInteraction.krw(satW)
                 / viscosityWBound;
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW)
                 / viscosityNwBound;
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
@@ -657,7 +664,7 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
                     (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (cellData.density(nPhaseIdx) + densityNwBound) : density_[nPhaseIdx];
         }
 
-        Scalar scalarPerm = permeability.two_norm();
+        const Scalar scalarPerm = permeability.two_norm();
 
         // calculate the gravity term
         Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
index 4056d8f703..c457759256 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh
@@ -42,6 +42,8 @@
 
 #include <dune/common/dynvector.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -96,7 +98,6 @@ class MimeticTwoPLocalStiffness: public LocalStiffness<TypeTag, 1>
     using CellData = GetPropType<TypeTag, Properties::CellData>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
 
 public:
     // define the number of components of your system, this is used outside
@@ -682,10 +683,13 @@ void MimeticTwoPLocalStiffness<TypeTag>::assembleElementMatrices(const Element&
                     PrimaryVariables boundValues(0.0);
                     problem_.dirichlet(boundValues, intersection);
 
-                    Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                            boundValues[saturationIdx]);
-                    Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
-                            boundValues[saturationIdx]);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    const Scalar krw = fluidMatrixInteraction.krw(boundValues[saturationIdx]);
+                    const Scalar krn = fluidMatrixInteraction.krn(boundValues[saturationIdx]);
 
                     switch (pressureType)
                     {
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
index 60585c6ccc..3de77fd91c 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh
@@ -44,6 +44,8 @@
 #include <dune/common/dynvector.hh>
 #include <dune/common/dynmatrix.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -98,7 +100,6 @@ class MimeticTwoPLocalStiffnessAdaptive: public LocalStiffness<TypeTag, 1>
     using CellData = GetPropType<TypeTag, Properties::CellData>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
 
     using IntersectionMapper = Dumux::IntersectionMapper<GridView>;
 
@@ -701,10 +702,13 @@ void MimeticTwoPLocalStiffnessAdaptive<TypeTag>::assembleElementMatrices(const E
                 PrimaryVariables boundValues(0.0);
                 problem_.dirichlet(boundValues, intersection);
 
-                Scalar krw = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                        boundValues[saturationIdx]);
-                Scalar krn = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
-                        boundValues[saturationIdx]);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                const Scalar krw = fluidMatrixInteraction.krw(boundValues[saturationIdx]);
+                const Scalar krn = fluidMatrixInteraction.krn(boundValues[saturationIdx]);
 
                 switch (pressureType)
                 {
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
index ca84c4a800..d4d10fbc0a 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/mimetic/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -57,7 +59,6 @@ class MimeticOperatorAssemblerTwoP: public CROperatorAssemblerTwoP<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
 
     using CellData = GetPropType<TypeTag, Properties::CellData>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -247,14 +248,18 @@ public:
                         PrimaryVariables boundValues(0.0);
                         problem.dirichlet(boundValues, intersection);
 
+                        // old material law interface is deprecated: Replace this by
+                        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                        // after the release of 3.3, when the deprecated interface is no longer supported
+                        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element);
+
                         if (velocityW[idxInInside] >= 0.)
                         {
                             mobilityW = cellData.mobility(wPhaseIdx);
                         }
                         else
                         {
-                            mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(element),
-                                    boundValues[saturationIdx]) / viscosityW;
+                            mobilityW = fluidMatrixInteraction.krw(boundValues[saturationIdx]) / viscosityW;
                         }
 
                         if (velocityNw[idxInInside] >= 0.)
@@ -263,8 +268,7 @@ public:
                         }
                         else
                         {
-                            mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(element),
-                                    boundValues[saturationIdx]) / viscosityNw;
+                            mobilityNw = fluidMatrixInteraction.krn(boundValues[saturationIdx]) / viscosityNw;
                         }
                     }
                     else
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
index d204e8af45..32135a2b63 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/mimetic/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -59,7 +61,6 @@ class MimeticOperatorAssemblerTwoPAdaptive : public CROperatorAssemblerTwoPAdapt
     using Element = typename GridView::template Codim<0>::Entity;
 
     using CellData = GetPropType<TypeTag, Properties::CellData>;
-    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -255,14 +256,18 @@ public:
                         PrimaryVariables boundValues(0.0);
                         problem.dirichlet(boundValues, intersection);
 
+                        // old material law interface is deprecated: Replace this by
+                        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                        // after the release of 3.3, when the deprecated interface is no longer supported
+                        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element);
+
                         if (velocityW[intersectionIdx] >= 0.)
                         {
                             mobilityW = cellData.mobility(wPhaseIdx);
                         }
                         else
                         {
-                            mobilityW = MaterialLaw::krw(problem.spatialParams().materialLawParams(element),
-                                boundValues[saturationIdx]) / viscosityW;
+                            mobilityW = fluidMatrixInteraction.krw(boundValues[saturationIdx]) / viscosityW;
                         }
 
                         if (velocityNw[intersectionIdx] >= 0.)
@@ -271,8 +276,7 @@ public:
                         }
                         else
                         {
-                            mobilityNw = MaterialLaw::krn(problem.spatialParams().materialLawParams(element),
-                                boundValues[saturationIdx]) / viscosityNw;
+                            mobilityNw = fluidMatrixInteraction.krn(boundValues[saturationIdx]) / viscosityNw;
                         }
                     }
                     else
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
index 92fa2c1259..9da440fb95 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressure.hh
@@ -32,6 +32,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operator.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimetic.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -67,7 +69,6 @@ template<class TypeTag> class MimeticPressure2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -515,13 +516,16 @@ void MimeticPressure2P<TypeTag>::updateMaterialLaws()
 
             CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-            Scalar satW = cellData.saturation(wPhaseIdx);
+            const Scalar satW = cellData.saturation(wPhaseIdx);
+
+            // old material law interface is deprecated: Replace this by
+            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+            // after the release of 3.3, when the deprecated interface is no longer supported
+            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
             // initialize mobilities
-            Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                    / viscosity_[wPhaseIdx];
-            Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                    / viscosity_[nPhaseIdx];
+            const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+            const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
             // initialize mobilities
             cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
index fb0f271f81..7ecf0b8839 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mimetic/pressureadaptive.hh
@@ -32,6 +32,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/operatoradaptive.hh>
 #include <dumux/porousmediumflow/2p/sequential/diffusion/mimetic/mimeticadaptive.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -67,7 +69,6 @@ template<class TypeTag> class MimeticPressure2PAdaptive
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -532,13 +533,16 @@ void MimeticPressure2PAdaptive<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
index 34c6c1fabb..d5491eaa97 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressure.hh
@@ -29,6 +29,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -82,7 +84,6 @@ class FvMpfaL2dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1664,19 +1665,20 @@ void FvMpfaL2dPressure2p<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -1840,15 +1842,18 @@ void FvMpfaL2dPressure2p<TypeTag>::updateMaterialLaws()
 
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
index e19b8d1740..6813decae3 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureadaptive.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -85,7 +87,6 @@ class FvMpfaL2dPressure2pAdaptive: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -2383,19 +2384,20 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -2556,17 +2558,20 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
index 59b4039824..15b59fadef 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocity.hh
@@ -27,6 +27,8 @@
 #include "2dpressure.hh"
 #include "2dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -61,7 +63,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2p: public FvMpfaL2dPress
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -453,8 +454,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -470,10 +477,8 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
index d4a7418692..87b3176ead 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurevelocityadaptive.hh
@@ -27,6 +27,8 @@
 #include "2dpressureadaptive.hh"
 #include "2dvelocityadaptive.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -63,7 +65,6 @@ template<class TypeTag> class FvMpfaL2dPressureVelocity2pAdaptive: public FvMpfa
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -589,8 +590,14 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -606,10 +613,8 @@ void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
index 269960826c..0e6df329cd 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocity.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
 #include "2dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -66,7 +68,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2p
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -786,21 +787,20 @@ void FvMpfaL2dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                             * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector < Scalar, numPhases
-                            > lambdaBound(
-                                    MaterialLaw::krw(
-                                            problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector <Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
index 2fddf014ca..85cbc3586b 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dvelocityadaptive.hh
@@ -63,7 +63,6 @@ template<class TypeTag> class FvMpfaL2dVelocity2pAdaptive : public FvMpfaL2dVelo
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
index 68a418328f..8eefe30623 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressure.hh
@@ -32,6 +32,8 @@
 #include "3dinteractionvolumecontainer.hh"
 #include "3dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -85,7 +87,6 @@ class FvMpfaL3dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -2351,18 +2352,20 @@ void FvMpfaL3dPressure2p<TypeTag>::assembleBoundaryInteractionVolume(Interaction
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                                                     problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                         * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector<Scalar, numPhases>
-                      lambdaBound(MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                                              problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -2489,14 +2492,19 @@ void FvMpfaL3dPressure2p<TypeTag>::updateMaterialLaws()
 
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar mobilityW = fluidMatrixInteraction.krw(satW)
             / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
+        Scalar mobilityNw = fluidMatrixInteraction.krn(satW)
             / viscosity_[nPhaseIdx];
 
         // initialize mobilities
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
index a5129ed870..ef2aef01af 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressureadaptive.hh
@@ -89,7 +89,6 @@ class FvMpfaL3dPressure2pAdaptive: public FvMpfaL3dPressure2p<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
index 025df282e2..595ff582d9 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocity.hh
@@ -28,6 +28,8 @@
 #include "3dpressure.hh"
 #include "3dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -77,7 +79,6 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2p: public FvMpfaL3dPress
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
     using Intersection = typename GridView::Intersection;
@@ -452,8 +453,14 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -469,10 +476,8 @@ void FvMpfaL3dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
index 5782b15acf..84a8ba579a 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dpressurevelocityadaptive.hh
@@ -27,6 +27,8 @@
 #include "3dpressureadaptive.hh"
 #include "3dvelocityadaptive.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -79,7 +81,6 @@ template<class TypeTag> class FvMpfaL3dPressureVelocity2pAdaptive: public FvMpfa
     using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
     using Intersection = typename GridView::Intersection;
@@ -553,8 +554,14 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -570,10 +577,8 @@ void FvMpfaL3dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(c
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-        / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-        / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
index 8aef1f4cb8..818f50b71b 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocity.hh
@@ -29,6 +29,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include "3dtransmissibilitycalculator.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -63,7 +65,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2p
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1972,19 +1973,20 @@ void FvMpfaL3dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                 }
 
-                Scalar pcBound = MaterialLaw::pc(
-                        problem_.spatialParams().materialLawParams(element), satWBound);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                         * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                 pcBound += gravityDiffBound;
 
-                Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                        MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                satWBound));
-                lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                        problem_.spatialParams().materialLawParams(element), satWBound);
+                Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
index 9aa9447370..1afe847c7f 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/3dvelocityadaptive.hh
@@ -65,7 +65,6 @@ template<class TypeTag> class FvMpfaL3dVelocity2pAdaptive: public FvMpfaL3dVeloc
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
index fe1b2e9970..97f3be5fe5 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressure.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -77,7 +79,6 @@ class FvMpfaO2dPressure2p: public FVPressure<TypeTag>
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -1763,19 +1764,20 @@ void FvMpfaO2dPressure2p<TypeTag>::assemble()
 
                             }
 
-                            Scalar pcBound = MaterialLaw::pc(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            // old material law interface is deprecated: Replace this by
+                            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                            // after the release of 3.3, when the deprecated interface is no longer supported
+                            const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                            Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                             Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                                     * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                             pcBound += gravityDiffBound;
 
-                            Dune::FieldVector<Scalar, numPhases> lambdaBound(
-                                    MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                            lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                                    problem_.spatialParams().materialLawParams(element), satWBound);
+                            Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                            lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                             lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                             lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
@@ -1938,17 +1940,20 @@ void FvMpfaO2dPressure2p<TypeTag>::updateMaterialLaws()
 
         CellData& cellData = problem_.variables().cellData(eIdxGlobal);
 
-        Scalar satW = cellData.saturation(wPhaseIdx);
+        const Scalar satW = cellData.saturation(wPhaseIdx);
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
index a38f5e39d1..eb00f5a077 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressurevelocity.hh
@@ -27,6 +27,8 @@
 #include "2dpressure.hh"
 #include "2dvelocity.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -60,7 +62,6 @@ template<class TypeTag> class FvMpfaO2dPressureVelocity2p: public FvMpfaO2dPress
     using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -454,8 +455,14 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             satW = cellData.saturation(wPhaseIdx);
         }
 
-        Scalar pressBound = boundValues[pressureIdx];
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        const Scalar pressBound = boundValues[pressureIdx];
+
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satW);
 
         //determine phase pressures from primary pressure variable
         Scalar pressWBound = 0;
@@ -471,10 +478,8 @@ void FvMpfaO2dPressureVelocity2p<TypeTag>::calculateVelocityOnBoundary(const Int
             pressNwBound = pressBound;
         }
 
-        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[wPhaseIdx];
-        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
-                / viscosity_[nPhaseIdx];
+        const Scalar lambdaWBound = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar lambdaNwBound = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
         Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
diff --git a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
index b9ede8717a..36a08b2bd7 100644
--- a/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
+++ b/dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dvelocity.hh
@@ -30,6 +30,8 @@
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/mpfa/ointeractionvolume.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -65,7 +67,6 @@ template<class TypeTag> class FvMpfaO2dVelocity2P
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -619,21 +620,20 @@ void FvMpfaO2dVelocity2P<TypeTag>::calculateBoundaryInteractionVolumeVelocity(In
 
                     }
 
-                    Scalar pcBound = MaterialLaw::pc(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    // old material law interface is deprecated: Replace this by
+                    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+                    // after the release of 3.3, when the deprecated interface is no longer supported
+                    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+                    Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
                     Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
                             * (density_[nPhaseIdx] - density_[wPhaseIdx]);
 
                     pcBound += gravityDiffBound;
 
-                    Dune::FieldVector < Scalar, numPhases
-                            > lambdaBound(
-                                    MaterialLaw::krw(
-                                            problem_.spatialParams().materialLawParams(element),
-                                            satWBound));
-                    lambdaBound[nPhaseIdx] = MaterialLaw::krn(
-                            problem_.spatialParams().materialLawParams(element), satWBound);
+                    Dune::FieldVector < Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
+                    lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
                     lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
                     lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
 
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
index 4e50446756..dfdc2b3ae0 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/capillarydiffusion.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/diffusivepart.hh>
 #include "properties.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -54,7 +56,6 @@ private:
       using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
       using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-      using MaterialLaw = typename SpatialParams::MaterialLaw;
 
       using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
       using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -112,6 +113,11 @@ public:
         Scalar mobilityWI = 0;
         Scalar mobilityNwI = 0;
 
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
         if (preComput_)
         {
             mobilityWI = CellDataI.mobility(wPhaseIdx);
@@ -123,9 +129,9 @@ public:
             fluidState.setPressure(wPhaseIdx, referencePressure);
             fluidState.setPressure(nPhaseIdx, referencePressure);
             fluidState.setTemperature(temperature);
-            mobilityWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satI);
+            mobilityWI = fluidMatrixInteraction.krw(satI);
             mobilityWI /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-            mobilityNwI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satI);
+            mobilityNwI = fluidMatrixInteraction.krn(satI);
             mobilityNwI /= FluidSystem::viscosity(fluidState, nPhaseIdx);
         }
 
@@ -172,9 +178,14 @@ public:
                 fluidState.setPressure(nPhaseIdx, referencePressure);
                 fluidState.setTemperature(temperature);
 
-                mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteractionNeighbor = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+                mobilityWJ = fluidMatrixInteractionNeighbor.krw(satJ);
                 mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-                mobilityNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                mobilityNwJ = fluidMatrixInteractionNeighbor.krn(satJ);
                 mobilityNwJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
             }
             Scalar mobilityWMean = 0.5*(mobilityWI + mobilityWJ);
@@ -207,9 +218,9 @@ public:
             fluidState.setPressure(wPhaseIdx, referencePressure);
             fluidState.setPressure(nPhaseIdx, referencePressure);
             fluidState.setTemperature(temperature);
-            mobilityWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satJ);
+            mobilityWJ = fluidMatrixInteraction.krw(satJ);
             mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
-            mobilityNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satJ);
+            mobilityNwJ = fluidMatrixInteraction.krn(satJ);
             mobilityNwJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
 
             Scalar mobWMean = 0.5 * (mobilityWI + mobilityWJ);
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
index 0f43ac4153..4750fe67ab 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh
@@ -28,6 +28,8 @@
 #include <dumux/porousmediumflow/sequential/impetproperties.hh>
 #include "evalcflflux.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -44,7 +46,6 @@ private:
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -357,7 +358,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
     Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
     Scalar lambdaNwI = cellDataI.mobility(nPhaseIdx);
 
-    Scalar dpc_dsI = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(element), satI);
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+    // after the release of 3.3, when the deprecated interface is no longer supported
+    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+    const Scalar dpc_dsI = fluidMatrixInteraction.dpc_dsw(satI);
 
     const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
 
@@ -407,7 +413,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
             Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
             Scalar lambdaNwJ = cellDataI.mobility(nPhaseIdx);
 
-            Scalar dpc_dsJ = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+            // old material law interface is deprecated: Replace this by
+            // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+            // after the release of 3.3, when the deprecated interface is no longer supported
+            const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+            const Scalar dpc_dsJ = fluidMatrixInteractionNeighbor.dpc_dsw(satJ);
 
             // compute vectorized permeabilities
             DimVector permeability(0);
@@ -448,8 +459,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), abs(satPlus)) / viscosityW;
-            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), abs(satMinus)) / viscosityW;
+            Scalar dLambdaWDs = fluidMatrixInteractionNeighbor.krw(abs(satPlus)) / viscosityW;
+            dLambdaWDs -= fluidMatrixInteractionNeighbor.krw(abs(satMinus)) / viscosityW;
             dLambdaWDs /= (ds);
 
             if (upwindNwI)
@@ -471,8 +482,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satPlus) / viscosityNw;
-            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satMinus) / viscosityNw;
+            Scalar dLambdaNwDs = fluidMatrixInteractionNeighbor.krn(satPlus) / viscosityNw;
+            dLambdaNwDs -= fluidMatrixInteractionNeighbor.krn(satMinus) / viscosityNw;
             dLambdaNwDs /= (ds);
 
             Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
@@ -559,13 +570,13 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                     case pw:
                     {
                         potWBound = bcValues[eqIdxPress] + density_[wPhaseIdx] * gdeltaZ;
-                        potNwBound = bcValues[eqIdxPress] + MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satWBound)
+                        potNwBound = bcValues[eqIdxPress] + fluidMatrixInteraction.pc(satWBound)
                                                           + density_[nPhaseIdx] * gdeltaZ;
                         break;
                     }
                     case pn:
                     {
-                        potWBound = bcValues[eqIdxPress] - MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                        potWBound = bcValues[eqIdxPress] - fluidMatrixInteraction.pc(satWBound)
                                                          + density_[wPhaseIdx] * gdeltaZ;
                         potNwBound = bcValues[eqIdxPress] + density_[nPhaseIdx] * gdeltaZ;
                         break;
@@ -602,12 +613,12 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
 
                 if (hasPotWBound && !hasPotNwBound)
                 {
-                    potNwBound = potWBound + MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                    potNwBound = potWBound + fluidMatrixInteraction.pc(satWBound)
                                            + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
                 }
                 else if (!hasPotWBound && hasPotNwBound)
                 {
-                   potWBound = potNwBound - MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),satWBound)
+                   potWBound = potNwBound - fluidMatrixInteraction.pc(satWBound)
                                           + (density_[nPhaseIdx] - density_[wPhaseIdx]) * gdeltaZ;
                 }
             }
@@ -651,7 +662,7 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 return;
             }
 
-            Scalar dpc_dsBound = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(element), satWBound);
+            const Scalar dpc_dsBound = fluidMatrixInteraction.dpc_dsw(satWBound);
 
             Scalar lambdaWBound = 0;
             Scalar lambdaNwBound = 0;
@@ -666,8 +677,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
             Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
             Scalar viscosityNwBound =
                 FluidSystem::viscosity(fluidState, nPhaseIdx);
-            lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satWBound) / viscosityWBound;
-            lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satWBound) / viscosityNwBound;
+            lambdaWBound = fluidMatrixInteraction.krw(satWBound) / viscosityWBound;
+            lambdaNwBound = fluidMatrixInteraction.krn(satWBound) / viscosityNwBound;
 
             Scalar satUpw = 0;
             using std::max;
@@ -690,8 +701,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satPlus) / viscosityW;
-            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satMinus) / viscosityW;
+            Scalar dLambdaWDs = fluidMatrixInteraction.krw(satPlus) / viscosityW;
+            dLambdaWDs -= fluidMatrixInteraction.krw(satMinus) / viscosityW;
             dLambdaWDs /= (ds);
 
             if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
@@ -713,8 +724,8 @@ void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw,
                 ds += epsDerivative_;
             }
 
-            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satPlus) / viscosityNw;
-            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satMinus) / viscosityNw;
+            Scalar dLambdaNwDs = fluidMatrixInteraction.krn(satPlus) / viscosityNw;
+            dLambdaNwDs -= fluidMatrixInteraction.krn(satMinus) / viscosityNw;
             dLambdaNwDs /= (ds);
 
             Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
index ae2e2ac4ec..d249715a25 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxdefault.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/sequential/impetproperties.hh>
 #include "evalcflflux.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -234,11 +236,16 @@ private:
 template<class TypeTag>
 typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCflFluxFunction(const Element& element)
 {
-    Scalar residualSatW = problem_.spatialParams().materialLawParams(element).swr();
-    Scalar residualSatNw = problem_.spatialParams().materialLawParams(element).snr();
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+    // after the release of 3.3, when the deprecated interface is no longer supported
+    const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+    const Scalar residualSatW = fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr();
+    const Scalar residualSatNw = fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr();
 
     // compute dt restriction
-    Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNw;
+    const Scalar volumeCorrectionFactor = 1 - residualSatW - residualSatNw;
     Scalar volumeCorrectionFactorOutW = 0;
     Scalar volumeCorrectionFactorOutNw = 0;
 
@@ -273,7 +280,7 @@ typename EvalCflFluxDefault<TypeTag>::Scalar EvalCflFluxDefault<TypeTag>::getCfl
     }
 
     //determine timestep
-    Scalar cFLFluxFunction = min(cFLFluxIn, cFLFluxOut);
+    const Scalar cFLFluxFunction = min(cFLFluxIn, cFLFluxOut);
 
     return cFLFluxFunction;
 }
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
index 59cb77ff6a..2e2eab7375 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/gravitypart.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/convectivepart.hh>
 #include "properties.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPModel
@@ -54,7 +56,6 @@ private:
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
@@ -100,6 +101,11 @@ public:
         Scalar lambdaWJ = 0;
         Scalar lambdaNwJ = 0;
 
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
         if (preComput_)
         {
             lambdaWI=cellDataI.mobility(wPhaseIdx);
@@ -107,9 +113,9 @@ public:
         }
         else
         {
-            lambdaWI = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satI);
+            lambdaWI = fluidMatrixInteraction.krw(satI);
             lambdaWI /= viscosity_[wPhaseIdx];
-            lambdaNwI = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satI);
+            lambdaNwI = fluidMatrixInteraction.krn(satI);
             lambdaNwI /= viscosity_[nPhaseIdx];
         }
 
@@ -144,9 +150,14 @@ public:
             }
             else
             {
-                lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                // old material law interface is deprecated: Replace this by
+                // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(neighbor.geometry().center());
+                // after the release of 3.3, when the deprecated interface is no longer supported
+                const auto fluidMatrixInteractionNeighbor = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), neighbor);
+
+                lambdaWJ = fluidMatrixInteractionNeighbor.krw(satJ);
                 lambdaWJ /= viscosity_[wPhaseIdx];
-                lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(neighbor), satJ);
+                lambdaNwJ = fluidMatrixInteractionNeighbor.krn(satJ);
                 lambdaNwJ /= viscosity_[nPhaseIdx];
             }
 
@@ -164,9 +175,9 @@ public:
             distVec = intersection.geometry().center() - element.geometry().center();
 
             //calculate lambda_n*f_w at the boundary
-            lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satJ);
+            lambdaWJ = fluidMatrixInteraction.krw(satJ);
             lambdaWJ /= viscosity_[wPhaseIdx];
-            lambdaNwJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satJ);
+            lambdaNwJ = fluidMatrixInteraction.krn(satJ);
             lambdaNwJ /= viscosity_[nPhaseIdx];
 
             //If potential is zero always take value from the boundary!
diff --git a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
index d8331db506..da29f8af87 100644
--- a/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
+++ b/dumux/porousmediumflow/2p/sequential/transport/cellcentered/saturation.hh
@@ -27,6 +27,8 @@
 #include <dumux/porousmediumflow/2p/sequential/transport/properties.hh>
 #include <dumux/porousmediumflow/sequential/cellcentered/transport.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 
 /*!
@@ -92,7 +94,6 @@ class FVSaturation2P: public FVTransport<TypeTag>
     using GravityFlux = GetPropType<TypeTag, Properties::GravityFlux>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
 
@@ -803,7 +804,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         }
         }
 
-        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(elementI), satWBound);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(elementI.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), elementI);
+
+        const Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
 
         Scalar lambdaW = 0;
         Scalar lambdaNw = 0;
@@ -821,12 +827,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         {
             if (compressibility_)
             {
-                lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaW = fluidMatrixInteraction.krw(satWBound)
                         / FluidSystem::viscosity(cellDataI.fluidState(), wPhaseIdx);
             }
             else
             {
-                lambdaW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaW = fluidMatrixInteraction.krw(satWBound)
                         / viscosity_[wPhaseIdx];
             }
         }
@@ -843,12 +849,12 @@ void FVSaturation2P<TypeTag>::getFluxOnBoundary(Scalar& update, const Intersecti
         {
             if (compressibility_)
             {
-                lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaNw = fluidMatrixInteraction.krn(satWBound)
                         / FluidSystem::viscosity(cellDataI.fluidState(), nPhaseIdx);
             }
             else
             {
-                lambdaNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(elementI), satWBound)
+                lambdaNw = fluidMatrixInteraction.krn(satWBound)
                         / viscosity_[nPhaseIdx];
             }
         }
@@ -1196,13 +1202,18 @@ void FVSaturation2P<TypeTag>::updateMaterialLaws()
         //determine phase saturations from primary saturation variable
         Scalar satW = cellData.saturation(wPhaseIdx);
 
-        Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(elementI.geometry().center());
+        // after the release of 3.3, when the deprecated interface is no longer supported
+        const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
+
+        const Scalar pc = fluidMatrixInteraction.pc(satW);
 
         cellData.setCapillaryPressure(pc);
 
         // initialize mobilities
-        Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[wPhaseIdx];
-        Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW) / viscosity_[nPhaseIdx];
+        const Scalar mobilityW = fluidMatrixInteraction.krw(satW) / viscosity_[wPhaseIdx];
+        const Scalar mobilityNw = fluidMatrixInteraction.krn(satW) / viscosity_[nPhaseIdx];
 
         // initialize mobilities
         cellData.setMobility(wPhaseIdx, mobilityW);
-- 
GitLab