From 85ef0f1d293cfc68e4651b68588b426de237f863 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 26 Oct 2020 11:23:03 +0100
Subject: [PATCH] [2p2c][sequential] Adapt headers

---
 .../2p2c/sequential/fv2dtransportadaptive.hh  | 11 ++++-
 .../2p2c/sequential/fv3dpressureadaptive.hh   |  1 -
 .../2p2c/sequential/fv3dtransportadaptive.hh  | 12 +++--
 .../2p2c/sequential/fvpressure.hh             | 30 +++++++-----
 .../sequential/fvpressurecompositional.hh     | 25 +++++-----
 .../2p2c/sequential/fvpressuremultiphysics.hh | 45 +++++++++--------
 .../2p2c/sequential/fvtransport.hh            | 49 ++++++++++++-------
 7 files changed, 105 insertions(+), 68 deletions(-)

diff --git a/dumux/porousmediumflow/2p2c/sequential/fv2dtransportadaptive.hh b/dumux/porousmediumflow/2p2c/sequential/fv2dtransportadaptive.hh
index b7937fb722..37754d5bde 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fv2dtransportadaptive.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fv2dtransportadaptive.hh
@@ -33,6 +33,8 @@
 #include "adaptiveproperties.hh"
 #include "fvtransport.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -366,13 +368,18 @@ void FV2dTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2
     Scalar pressI = problem().pressureModel().pressure(globalIdxI);
     Scalar pcI = cellDataI.capillaryPressure();
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     PhaseVector SmobI(0.);
     using std::max;
     SmobI[wPhaseIdx] = max((cellDataI.saturation(wPhaseIdx)
-                            - problem().spatialParams().materialLawParams(elementI).swr())
+                            - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
                             , 1e-2);
     SmobI[nPhaseIdx] = max((cellDataI.saturation(nPhaseIdx)
-                                - problem().spatialParams().materialLawParams(elementI).snr())
+                                - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
                             , 1e-2);
 
     Scalar densityWI (0.), densityNWI(0.);
diff --git a/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh b/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
index e4cc63dfa1..2092c28fda 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fv3dpressureadaptive.hh
@@ -90,7 +90,6 @@ template<class TypeTag> class FV3dPressure2P2CAdaptive
     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;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
diff --git a/dumux/porousmediumflow/2p2c/sequential/fv3dtransportadaptive.hh b/dumux/porousmediumflow/2p2c/sequential/fv3dtransportadaptive.hh
index 4e505d0e70..e18141fe40 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fv3dtransportadaptive.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fv3dtransportadaptive.hh
@@ -33,6 +33,8 @@
 #include "adaptiveproperties.hh"
 #include "fvtransport.hh"
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -57,7 +59,6 @@ class FV3dTransport2P2CAdaptive : public FVTransport2P2C<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;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -364,13 +365,18 @@ void FV3dTransport2P2CAdaptive<TypeTag>::getMpfaFlux(Dune::FieldVector<Scalar, 2
     Scalar pressI = problem().pressureModel().pressure(globalIdxI);
     Scalar pcI = cellDataI.capillaryPressure();
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     PhaseVector SmobI(0.);
     using std::max;
     SmobI[wPhaseIdx] = max((cellDataI.saturation(wPhaseIdx)
-                            - problem().spatialParams().materialLawParams(elementI).swr())
+                            - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
                             , 1e-2);
     SmobI[nPhaseIdx] = max((cellDataI.saturation(nPhaseIdx)
-                                - problem().spatialParams().materialLawParams(elementI).snr())
+                                - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
                             , 1e-2);
 
     Scalar densityWI (0.), densityNWI(0.);
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
index 0b14ff3425..ea035883dc 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressure.hh
@@ -38,6 +38,8 @@
 #include <dumux/io/vtkmultiwriter.hh>
 #include <dumux/porousmediumflow/2p2c/sequential/properties.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -79,7 +81,6 @@ template<class TypeTag> class FVPressure2P2C
     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;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -677,6 +678,11 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
     PhaseVector pressBC(0.);
     Scalar pcBound (0.);
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     /**********         Dirichlet Boundary        *************/
     if (bcType.isDirichlet(Indices::pressureEqIdx))
     {
@@ -754,11 +760,9 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
             else if(getPropValue<TypeTag, Properties::BoundaryMobility>() == Indices::permDependent)
             {
                 lambdaWBound
-                    = MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI),
-                            BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
+                    = fluidMatrixInteraction.krw(BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
                 lambdaNWBound
-                    = MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI),
-                            BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
+                    = fluidMatrixInteraction.krn(BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
             }
             // get average density
             Scalar rhoMeanW = 0.5 * (cellDataI.density(wPhaseIdx) + densityWBound);
@@ -953,6 +957,13 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
 
     PhaseVector pressure;
     CompositionalFlash<Scalar, FluidSystem> flashSolver;
+
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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(getPropValue<TypeTag, Properties::EnableCapillarity>()) // iterate capillary pressure and saturation
     {
         unsigned int maxiter = 6;
@@ -981,8 +992,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
 
             // calculate new pc
             Scalar oldPc = pc;
-            pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
-                                 fluidState.saturation(wPhaseIdx));
+            pc = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
 
             if (fabs(oldPc-pc)<10 && iter != 0)
                 break;
@@ -999,11 +1009,9 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
     }
 
     // initialize mobilities
-    cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem().spatialParams().materialLawParams(element),
-                                                     fluidState.saturation(wPhaseIdx))
+    cellData.setMobility(wPhaseIdx, fluidMatrixInteraction.krw(fluidState.saturation(wPhaseIdx))
                 / cellData.viscosity(wPhaseIdx));
-    cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem().spatialParams().materialLawParams(element),
-                                                     fluidState.saturation(wPhaseIdx))
+    cellData.setMobility(nPhaseIdx, fluidMatrixInteraction.krn(fluidState.saturation(wPhaseIdx))
                 / cellData.viscosity(nPhaseIdx));
 
     // determine volume mismatch between actual fluid volume and pore volume
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
index f99bd8dad9..e0734b718f 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressurecompositional.hh
@@ -34,6 +34,8 @@
 #include <dumux/porousmediumflow/2p2c/sequential/properties.hh>
 #include <dumux/io/vtkmultiwriter.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -78,9 +80,6 @@ template<class TypeTag> class FVPressureCompositional
 
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
-    ///@cond false
-    using MaterialLaw = typename GetPropType<TypeTag, Properties::SpatialParams>::MaterialLaw;
-    ///@endcond
 
     using CellData = GetPropType<TypeTag, Properties::CellData>;
     enum
@@ -547,6 +546,11 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
         FluidState& fluidState = cellData.manipulateFluidState();
         CompositionalFlash<Scalar, FluidSystem> flashSolver;
 
+        // 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);
+
         // initial conditions
         PhaseVector pressure(0.);
         Scalar sat_0=0.;
@@ -579,8 +583,7 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
                 Scalar pc=0.;
                 if(getPropValue<TypeTag, Properties::EnableCapillarity>())
                 {
-                    pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),
-                                    sat_0);
+                    pc = fluidMatrixInteraction.pc(sat_0);
                 }
                 else
                     pc = 0.;
@@ -640,16 +643,14 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
                         Scalar oldPc = pc;
                         //update with better pressures
                         flashSolver.concentrationFlash2p2c(fluidState, Z0, pressure, problem_.temperatureAtPos(globalPos));
-                        pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),
-                                            fluidState.saturation(wPhaseIdx));
+                        pc = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
                         // TODO: get right criterion, do output for evaluation
                         //converge criterion
                         using std::abs;
                         if (abs(oldPc - pc) < 10.0)
                             iter = maxiter;
 
-                        pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element),
-                                fluidState.saturation(wPhaseIdx));
+                        pc = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx));
                     }
                 }
                 else  // capillary pressure neglected
@@ -667,11 +668,9 @@ void FVPressureCompositional<TypeTag>::initialMaterialLaws(bool compositional)
         problem_.transportModel().totalConcentration(nCompIdx,eIdxGlobal) = cellData.massConcentration(nCompIdx);
 
         // initialize mobilities
-        cellData.setMobility(wPhaseIdx, MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),
-                                                         fluidState.saturation(wPhaseIdx))
+        cellData.setMobility(wPhaseIdx, fluidMatrixInteraction.krw(fluidState.saturation(wPhaseIdx))
                     / cellData.viscosity(wPhaseIdx));
-        cellData.setMobility(nPhaseIdx, MaterialLaw::krn(problem_.spatialParams().materialLawParams(element),
-                                                         fluidState.saturation(wPhaseIdx))
+        cellData.setMobility(nPhaseIdx, fluidMatrixInteraction.krn(fluidState.saturation(wPhaseIdx))
                     / cellData.viscosity(nPhaseIdx));
 
         // calculate perimeter used as weighting factor
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh b/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
index ed873d286d..2b6540239f 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvpressuremultiphysics.hh
@@ -31,6 +31,8 @@
 #include <dumux/parallel/vectorcommdatahandle.hh>
 #include <dumux/material/constraintsolvers/compositionalflash.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -75,7 +77,6 @@ class FVPressure2P2CMultiPhysics : public FVPressure2P2C<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;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -700,23 +701,23 @@ void FVPressure2P2CMultiPhysics<TypeTag>::get1pFluxOnBoundary(Dune::FieldVector<
                         Scalar lambdaBound = 0.;
                         switch (getPropValue<TypeTag, Properties::BoundaryMobility>())
                         {
-                        case Indices::satDependent:
+                            case Indices::satDependent:
                             {
-                                lambdaBound = BCfluidState.saturation(phaseIdx)
-                                    / viscosityBound;
-                            break;
+                                lambdaBound = BCfluidState.saturation(phaseIdx) / viscosityBound;
+                                break;
                             }
-                        case Indices::permDependent:
+                            case Indices::permDependent:
                             {
-                            if (phaseIdx == wPhaseIdx)
-                                lambdaBound = MaterialLaw::krw(
-                                    problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
-                                    / viscosityBound;
-                            else
-                                lambdaBound = MaterialLaw::krn(
-                                    problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
-                                    / viscosityBound;
-                            break;
+                                // old material law interface is deprecated: Replace this by
+                                // const auto& fluidMatrixInteraction = problem().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);
+
+                                if (phaseIdx == wPhaseIdx)
+                                    lambdaBound = fluidMatrixInteraction.krw(BCfluidState.saturation(wPhaseIdx)) / viscosityBound;
+                                else
+                                    lambdaBound = fluidMatrixInteraction.krn(BCfluidState.saturation(wPhaseIdx)) / viscosityBound;
+                                break;
                             }
                         }
                         Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + densityBound);
@@ -931,14 +932,18 @@ void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const El
     // acess the simple fluid state and prepare for manipulation
     auto& pseudoFluidState = cellData.manipulateSimpleFluidState();
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     // prepare phase pressure for fluid state
     // both phase pressures are necessary for the case 1p domain is assigned for
     // the next 2p subdomain
     PhaseVector pressure(0.);
     Scalar pc = 0;
     if(getPropValue<TypeTag, Properties::EnableCapillarity>())
-        pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(elementI),
-            ((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign sw = 1 if wPhase present, else 0
+        pc = fluidMatrixInteraction.pc(((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign sw = 1 if wPhase present, else 0
     if(pressureType == wPhaseIdx)
     {
         pressure[wPhaseIdx] = this->pressure(eIdxGlobal);
@@ -970,15 +975,13 @@ void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const El
     if(presentPhaseIdx == wPhaseIdx)
     {
         cellData.setMobility(wPhaseIdx,
-            MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
-                / cellData.viscosity(wPhaseIdx));
+            fluidMatrixInteraction.krw(pseudoFluidState.saturation(wPhaseIdx)) / cellData.viscosity(wPhaseIdx));
         cellData.setMobility(nPhaseIdx, 0.);
     }
     else
     {
         cellData.setMobility(nPhaseIdx,
-            MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
-                / cellData.viscosity(nPhaseIdx));
+            fluidMatrixInteraction.krn(pseudoFluidState.saturation(wPhaseIdx)) / cellData.viscosity(nPhaseIdx));
         cellData.setMobility(wPhaseIdx, 0.);
     }
 
diff --git a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
index 3bd9aa0be4..eff7ac11f3 100644
--- a/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/fvtransport.hh
@@ -35,6 +35,8 @@
 #include <dumux/common/math.hh>
 #include <dumux/parallel/vectorcommdatahandle.hh>
 
+#include <dumux/common/deprecated.hh>
+
 namespace Dumux {
 /*!
  * \ingroup SequentialTwoPTwoCModel
@@ -64,7 +66,6 @@ class FVTransport2P2C
     using Implementation = GetPropType<TypeTag, Properties::TransportModel>;
 
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
 
     using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
     using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
@@ -598,13 +599,18 @@ void FVTransport2P2C<TypeTag>::getFlux(ComponentVector& fluxEntries,
     Scalar pcI = cellDataI.capillaryPressure();
     DimMatrix K_I(problem().spatialParams().intrinsicPermeability(elementI));
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     PhaseVector SmobI(0.);
     using std::max;
     SmobI[wPhaseIdx] = max((cellDataI.saturation(wPhaseIdx)
-                            - problem().spatialParams().materialLawParams(elementI).swr())
+                            - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
                             , 1e-2);
     SmobI[nPhaseIdx] = max((cellDataI.saturation(nPhaseIdx)
-                                - problem().spatialParams().materialLawParams(elementI).snr())
+                                - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
                             , 1e-2);
 
     Scalar densityWI (0.), densityNWI(0.);
@@ -909,11 +915,16 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(ComponentVector& fluxEntries,
             K_I[axis][axis] = minimalBoundaryPermeability;
     }
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     Scalar SwmobI = max((cellDataI.saturation(wPhaseIdx)
-                            - problem().spatialParams().materialLawParams(elementI).swr())
+                            - fluidMatrixInteraction.pcSwCurve().effToAbsParams().swr())
                             , 1e-2);
     Scalar SnmobI = max((cellDataI.saturation(nPhaseIdx)
-                                - problem().spatialParams().materialLawParams(elementI).snr())
+                                - fluidMatrixInteraction.pcSwCurve().effToAbsParams().snr())
                             , 1e-2);
 
     Scalar densityWI (0.), densityNWI(0.);
@@ -1001,6 +1012,11 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(ComponentVector& fluxEntries,
         potential[wPhaseIdx] *= fabs(K * unitOuterNormal);
         potential[nPhaseIdx] *= fabs(K * unitOuterNormal);
 
+        // old material law interface is deprecated: Replace this by
+        // const auto& fluidMatrixInteraction = problem().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);
+
         // do upwinding for lambdas
         PhaseVector lambda(0.);
         if (potential[wPhaseIdx] >= 0.)
@@ -1010,9 +1026,7 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(ComponentVector& fluxEntries,
             if(getPropValue<TypeTag, Properties::BoundaryMobility>()==Indices::satDependent)
                 lambda[wPhaseIdx] = BCfluidState.saturation(wPhaseIdx) / viscosityWBound;
             else
-                lambda[wPhaseIdx] = MaterialLaw::krw(
-                        problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
-                        / viscosityWBound;
+                lambda[wPhaseIdx] = fluidMatrixInteraction.krw(BCfluidState.saturation(wPhaseIdx)) / viscosityWBound;
             }
         if (potential[nPhaseIdx] >= 0.)
             lambda[nPhaseIdx] = cellDataI.mobility(nPhaseIdx);
@@ -1021,9 +1035,7 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(ComponentVector& fluxEntries,
             if(getPropValue<TypeTag, Properties::BoundaryMobility>()==Indices::satDependent)
                 lambda[nPhaseIdx] = BCfluidState.saturation(nPhaseIdx) / viscosityNWBound;
             else
-                lambda[nPhaseIdx] = MaterialLaw::krn(
-                        problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
-                        / viscosityNWBound;
+                lambda[nPhaseIdx] = fluidMatrixInteraction.krn(BCfluidState.saturation(wPhaseIdx)) / viscosityNWBound;
             }
         // calculate and standardized velocity
 
@@ -1117,16 +1129,21 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
     PrimaryVariables primaryVariablesOnBoundary(0.);
     problem().dirichlet(primaryVariablesOnBoundary, intersection);
 
+    // old material law interface is deprecated: Replace this by
+    // const auto& fluidMatrixInteraction = problem().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);
+
     // read boundary type
     typename Indices::BoundaryFormulation bcType;
     problem().boundaryFormulation(bcType, intersection);
     if (bcType == Indices::saturation)
     {
         Scalar satBound = primaryVariablesOnBoundary[contiWEqIdx];
+
         if(getPropValue<TypeTag, Properties::EnableCapillarity>())
         {
-            Scalar pcBound = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
-                    satBound);
+            Scalar pcBound = fluidMatrixInteraction.pc(satBound);
             switch (pressureType)
             {
             case pw:
@@ -1159,8 +1176,7 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
 
         if(getPropValue<TypeTag, Properties::EnableCapillarity>())
         {
-            Scalar pcBound = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
-                    BCfluidState.saturation(wPhaseIdx));
+            Scalar pcBound = fluidMatrixInteraction.pc(BCfluidState.saturation(wPhaseIdx));
             int maxiter = 3;
             //start iteration loop
             for(int iter=0; iter < maxiter; iter++)
@@ -1189,8 +1205,7 @@ void FVTransport2P2C<TypeTag>::evalBoundary(GlobalPosition globalPosFace,
                 //update with better pressures
                 flashSolver.concentrationFlash2p2c(BCfluidState, Z0Bound, pressBound,
                         problem().temperatureAtPos(globalPosFace));
-                pcBound = MaterialLaw::pc(problem().spatialParams().materialLawParams(element),
-                        BCfluidState.saturation(wPhaseIdx));
+                pcBound = fluidMatrixInteraction.pc(BCfluidState.saturation(wPhaseIdx));
                 // TODO: get right criterion, do output for evaluation
                 //converge criterion
                 using std::abs;
-- 
GitLab