From 544c82dec427f93a432bb7fe2167cc2d16332ece Mon Sep 17 00:00:00 2001
From: Yonghui <huangyh56@gmail.com>
Date: Tue, 9 May 2017 00:36:16 +0200
Subject: [PATCH] Extract the two phase flow material property model to avoid
 code duplication

---
 .../CreateTwoPhaseFlowMaterialProperties.cpp  | 118 +++++++++++------
 .../CreateTwoPhaseFlowMaterialProperties.h    |  14 +-
 .../TwoPhaseFlowWithPPMaterialProperties.cpp  | 123 +++++++++++++-----
 .../TwoPhaseFlowWithPPMaterialProperties.h    | 103 ++++++++++-----
 4 files changed, 248 insertions(+), 110 deletions(-)

diff --git a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
index 21f9b1013ce..9d5849c5a5e 100644
--- a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
+++ b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
@@ -8,91 +8,135 @@
  */
 
 #include "CreateTwoPhaseFlowMaterialProperties.h"
+
 #include <logog/include/logog.hpp>
+
 #include "BaseLib/reorderVector.h"
 #include "MaterialLib/Fluid/FluidProperty.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/PropertyVector.h"
 #include "ProcessLib/Parameter/Parameter.h"
 #include "ProcessLib/Parameter/SpatialPosition.h"
+
 #include "TwoPhaseFlowWithPPMaterialProperties.h"
 
 namespace MaterialLib
 {
 namespace TwoPhaseFlowWithPP
 {
-std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
-CreateTwoPhaseFlowMaterialProperties(
+std::tuple<std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>,
+           BaseLib::ConfigTree>
+createTwoPhaseFlowMaterialProperties(
     BaseLib::ConfigTree const& config,
-    bool const has_material_ids,
-    MeshLib::PropertyVector<int> const& material_ids)
+    boost::optional<MeshLib::PropertyVector<int> const&>
+        material_ids)
 {
     DBUG("Reading material properties of two-phase flow process.");
 
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid}
-    auto const& fluid_config = config.getConfigSubtree("fluid");
+    //! \ogs_file_param{material__twophase_flow__material_property__fluid}
+    auto fluid_config = config.getConfigSubtree("fluid");
 
     // Get fluid properties
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid__liquid_density}
+    //! \ogs_file_param{material__twophase_flow__material_property__fluid__liquid_density}
     auto const& rho_conf = fluid_config.getConfigSubtree("liquid_density");
-    auto _liquid_density =
-        MaterialLib::Fluid::createFluidDensityModel(rho_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid__gas_density}
+    auto liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
+    //! \ogs_file_param{material__twophase_flow__material_property__fluid__gas_density}
     auto const& rho_gas_conf = fluid_config.getConfigSubtree("gas_density");
-    auto _gas_density =
+    auto gas_density =
         MaterialLib::Fluid::createFluidDensityModel(rho_gas_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid__liquid_viscosity}
+    //! \ogs_file_param{material__twophase_flow__material_property__fluid__liquid_viscosity}
     auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
-    auto _viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid__gas_viscosity}
+    auto liquid_viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+    //! \ogs_file_param{material__twophase_flow__material_property__fluid__gas_viscosity}
     auto const& mu_gas_conf = fluid_config.getConfigSubtree("gas_viscosity");
-    auto _gas_viscosity = MaterialLib::Fluid::createViscosityModel(mu_gas_conf);
+    auto gas_viscosity = MaterialLib::Fluid::createViscosityModel(mu_gas_conf);
 
     // Get porous properties
     std::vector<int> mat_ids;
-    std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
+    std::vector<int> mat_krel_ids;
+    std::vector<Eigen::MatrixXd> intrinsic_permeability_models;
     std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
-        _porosity_models;
+        porosity_models;
     std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
-        _storage_models;
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium}
+        storage_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+        capillary_pressure_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+        relative_permeability_models;
+
+    //! \ogs_file_param{material__twophase_flow__material_property__porous_medium}
     auto const& poro_config = config.getConfigSubtree("porous_medium");
-    //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium}
+    //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium}
     for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
     {
-        //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__id}
+        //! \ogs_file_attr{material__twophase_flow__material_property__porous_medium__porous_medium__id}
         auto const id = conf.getConfigAttributeOptional<int>("id");
         mat_ids.push_back(*id);
 
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__permeability}
+        //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__permeability}
         auto const& permeability_conf = conf.getConfigSubtree("permeability");
-        _intrinsic_permeability_models.emplace_back(
-            MaterialLib::PorousMedium::createPermeabilityModel(permeability_conf));
+        intrinsic_permeability_models.emplace_back(
+            MaterialLib::PorousMedium::createPermeabilityModel(
+                permeability_conf));
 
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__porosity}
+        //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__porosity}
         auto const& porosity_conf = conf.getConfigSubtree("porosity");
         auto n = MaterialLib::PorousMedium::createPorosityModel(porosity_conf);
-        _porosity_models.emplace_back(std::move(n));
+        porosity_models.emplace_back(std::move(n));
 
-        //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__storage}
+        //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__storage}
         auto const& storage_conf = conf.getConfigSubtree("storage");
         auto beta = MaterialLib::PorousMedium::createStorageModel(storage_conf);
-        _storage_models.emplace_back(std::move(beta));
+        storage_models.emplace_back(std::move(beta));
+
+        auto const& capillary_pressure_conf =
+            //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__capillary_pressure}
+            conf.getConfigSubtree("capillary_pressure");
+        auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel(
+            capillary_pressure_conf);
+        capillary_pressure_models.emplace_back(std::move(pc));
+
+        auto const& krel_config =
+            //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__relative_permeability}
+            conf.getConfigSubtree("relative_permeability");
+        for (
+            auto const& krel_conf :
+            //! \ogs_file_param{material__twophase_flow__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability}
+            krel_config.getConfigSubtreeList("relative_permeability"))
+        {
+            auto const krel_id =
+                //! \ogs_file_attr{material__twophase_flow__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability__id}
+                krel_conf.getConfigAttributeOptional<int>("id");
+            mat_krel_ids.push_back(*krel_id);
+            auto krel_n =
+                MaterialLib::PorousMedium::createRelativePermeabilityModel(
+                    krel_conf);
+            relative_permeability_models.emplace_back(std::move(krel_n));
+        }
+        BaseLib::reorderVector(relative_permeability_models, mat_krel_ids);
     }
 
-    BaseLib::reorderVector(_intrinsic_permeability_models, mat_ids);
-    BaseLib::reorderVector(_porosity_models, mat_ids);
-    BaseLib::reorderVector(_storage_models, mat_ids);
+    BaseLib::reorderVector(intrinsic_permeability_models, mat_ids);
+    BaseLib::reorderVector(porosity_models, mat_ids);
+    BaseLib::reorderVector(storage_models, mat_ids);
 
-    return std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>{
-        new TwoPhaseFlowWithPPMaterialProperties{
-            has_material_ids, material_ids, std::move(_liquid_density),
-            std::move(_viscosity), std::move(_gas_density),
-            std::move(_gas_viscosity), _intrinsic_permeability_models,
-            std::move(_porosity_models), std::move(_storage_models)}};
+    return std::forward_as_tuple(
+        std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>{
+            new TwoPhaseFlowWithPPMaterialProperties{
+                material_ids, std::move(liquid_density),
+                std::move(liquid_viscosity), std::move(gas_density),
+                std::move(gas_viscosity), std::move(intrinsic_permeability_models),
+                std::move(porosity_models), std::move(storage_models),
+                std::move(capillary_pressure_models),
+                std::move(relative_permeability_models)}},
+        std::move(fluid_config));
 }
 
 }  // end namespace
diff --git a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h
index bbf4f4309a9..19ebfd81205 100644
--- a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h
+++ b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h
@@ -10,11 +10,14 @@
 #pragma once
 
 #include <memory>
+#include <tuple>
+
 #include "MaterialLib/Fluid/FluidPropertyHeaders.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
-#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
+
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
 namespace BaseLib
 {
 class ConfigTree;
@@ -24,11 +27,12 @@ namespace MaterialLib
 {
 namespace TwoPhaseFlowWithPP
 {
-std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
-CreateTwoPhaseFlowMaterialProperties(
+std::tuple<std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>,
+           BaseLib::ConfigTree>
+createTwoPhaseFlowMaterialProperties(
     BaseLib::ConfigTree const& config,
-    bool const has_material_ids,
-    MeshLib::PropertyVector<int> const& material_ids);
+    boost::optional<MeshLib::PropertyVector<int> const&>
+        material_ids);
 
 }  // end namespace
 }  // end namespace
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
index 0b4c981f926..68fa49d7410 100644
--- a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
@@ -8,12 +8,17 @@
  */
 
 #include "TwoPhaseFlowWithPPMaterialProperties.h"
+
+#include <boost/math/special_functions/pow.hpp>
 #include <logog/include/logog.hpp>
 #include <utility>
+
 #include "BaseLib/reorderVector.h"
 #include "MaterialLib/Fluid/FluidProperty.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/PropertyVector.h"
@@ -25,46 +30,46 @@ namespace MaterialLib
 namespace TwoPhaseFlowWithPP
 {
 TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
-    bool const has_material_ids,
-    MeshLib::PropertyVector<int> const& material_ids,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        liquid_density,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        viscosity,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        gas_density,
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-        gas_viscosity,
-    std::vector<Eigen::MatrixXd>
-        intrinsic_permeability_models,
+    boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& liquid_density,
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& liquid_viscosity,
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& gas_density,
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& gas_viscosity,
+    std::vector<Eigen::MatrixXd>&& intrinsic_permeability_models,
     std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
         porosity_models,
     std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
-        storage_models)
-    : _has_material_ids(has_material_ids),
+        storage_models,
+    std::vector<std::unique_ptr<
+        MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+        capillary_pressure_models,
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+        relative_permeability_models)
+    : _material_ids(material_ids),
       _liquid_density(std::move(liquid_density)),
-      _viscosity(std::move(viscosity)),
+      _liquid_viscosity(std::move(liquid_viscosity)),
       _gas_density(std::move(gas_density)),
       _gas_viscosity(std::move(gas_viscosity)),
-      _material_ids(material_ids),
       _intrinsic_permeability_models(std::move(intrinsic_permeability_models)),
       _porosity_models(std::move(porosity_models)),
-      _storage_models(std::move(storage_models))
+      _storage_models(std::move(storage_models)),
+      _capillary_pressure_models(std::move(capillary_pressure_models)),
+      _relative_permeability_models(std::move(relative_permeability_models))
 {
     DBUG("Create material properties for Two-Phase flow with PP model.");
 }
 
-void TwoPhaseFlowWithPPMaterialProperties::setMaterialID(
-    const ProcessLib::SpatialPosition& pos)
+int TwoPhaseFlowWithPPMaterialProperties::getMaterialID(
+    const std::size_t element_id) const
 {
-    if (!_has_material_ids)
+    if (!_material_ids)
     {
-        _current_material_id = 0;
-        return;
+        return 0;
     }
 
-    assert(pos.getElementID().get() < _material_ids.size());
-    _current_material_id = _material_ids[pos.getElementID().get()];
+    assert(element_id < _material_ids->size());
+    return (*_material_ids)[element_id];
 }
 
 double TwoPhaseFlowWithPPMaterialProperties::getLiquidDensity(
@@ -85,15 +90,15 @@ double TwoPhaseFlowWithPPMaterialProperties::getGasDensity(const double p,
     return _gas_density->getValue(vars);
 }
 
-double TwoPhaseFlowWithPPMaterialProperties::getDerivGasDensity(
+double TwoPhaseFlowWithPPMaterialProperties::getGasDensityDerivative(
     const double p, const double T) const
 {
     ArrayType vars;
     vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
     vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
 
-    return _gas_density->getdValue(
-        vars, MaterialLib::Fluid::PropertyVariableType::p);
+    return _gas_density->getdValue(vars,
+                                   MaterialLib::Fluid::PropertyVariableType::p);
 }
 double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity(
     const double p, const double T) const
@@ -101,7 +106,7 @@ double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity(
     ArrayType vars;
     vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
     vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
-    return _viscosity->getValue(vars);
+    return _liquid_viscosity->getValue(vars);
 }
 
 double TwoPhaseFlowWithPPMaterialProperties::getGasViscosity(
@@ -114,21 +119,67 @@ double TwoPhaseFlowWithPPMaterialProperties::getGasViscosity(
 }
 
 Eigen::MatrixXd const& TwoPhaseFlowWithPPMaterialProperties::getPermeability(
-    const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
-    const int /*dim*/) const
+    const int material_id, const double /*t*/,
+    const ProcessLib::SpatialPosition& /*pos*/, const int /*dim*/) const
 {
-    return _intrinsic_permeability_models[_current_material_id];
+    return _intrinsic_permeability_models[material_id];
 }
 
 double TwoPhaseFlowWithPPMaterialProperties::getPorosity(
-    const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
-    const double /*p*/, const double T, const double porosity_variable) const
+    const int material_id, const double /*t*/,
+    const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+    const double T, const double porosity_variable) const
+{
+    return _porosity_models[material_id]->getValue(porosity_variable, T);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getSaturation(
+    const int material_id, const double /*t*/,
+    const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+    const double /*T*/, const double pc) const
 {
-    const double porosity =
-        _porosity_models[_current_material_id]->getValue(porosity_variable, T);
+    return _capillary_pressure_models[material_id]->getSaturation(pc);
+}
 
-    return porosity;
+double TwoPhaseFlowWithPPMaterialProperties::getCapillaryPressure(
+    const int material_id, const double /*t*/,
+    const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+    const double /*T*/, const double saturation) const
+{
+    return _capillary_pressure_models[material_id]->getCapillaryPressure(
+        saturation);
 }
 
+double TwoPhaseFlowWithPPMaterialProperties::getSaturationDerivative(
+    const int material_id, const double /*t*/,
+    const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+    const double /*T*/, const double saturation) const
+{
+    const double dpcdsw =
+        _capillary_pressure_models[material_id]->getdPcdS(saturation);
+    return 1 / dpcdsw;
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getNonwetRelativePermeability(
+    const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+    const double /*p*/, const double /*T*/, const double saturation) const
+{
+    if (saturation < 0.)
+        return 1.0;
+    else if (saturation > 1)
+        return 0.0;
+    return boost::math::pow<3>(1 - saturation);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getWetRelativePermeability(
+    const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+    const double /*p*/, const double /*T*/, const double saturation) const
+{
+    if (saturation < 0)
+        return 0.0;
+    else if (saturation > 1)
+        return 1.0;
+    return boost::math::pow<3>(saturation);
+}
 }  // end of namespace
 }  // end of namespace
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
index 977cc26e61c..83f7ed48358 100644
--- a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
@@ -12,10 +12,15 @@
 #include <iostream>
 #include <memory>
 #include <vector>
+
 #include "MaterialLib/Fluid/FluidPropertyHeaders.h"
 #include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 #include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
 #include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 
 // TODO
@@ -36,67 +41,101 @@ namespace MaterialLib
 {
 namespace TwoPhaseFlowWithPP
 {
+/** This class has a collection of material properties for two-phase flow with
+* PP model
+*  and it makes description of the material properties for two-phase condition,
+*  i.e. the gas/liquid density and viscosity models, respectively,
+*  the relative permeability models with respect to two phases,
+*  the capillary pressure-saturation relationships.
+*  It generally provides the computation of the PDE coefficients for two-phase
+* flow.
+*/
 class TwoPhaseFlowWithPPMaterialProperties
 {
 public:
     using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
 
     TwoPhaseFlowWithPPMaterialProperties(
-        bool const has_material_ids,
-        MeshLib::PropertyVector<int> const& material_ids,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+        boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
             liquid_density,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
-            viscosity,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
+            liquid_viscosity,
+        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
             gas_density,
-        std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+        std::unique_ptr<MaterialLib::Fluid::FluidProperty>&&
             gas_viscosity,
-        std::vector<Eigen::MatrixXd>
+        std::vector<Eigen::MatrixXd>&&
             intrinsic_permeability_models,
         std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
             porosity_models,
         std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
-            storage_models);
+            storage_models,
+        std::vector<std::unique_ptr<
+            MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+            capillary_pressure_models,
+        std::vector<
+            std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+            relative_permeability_models);
 
-    void setMaterialID(const ProcessLib::SpatialPosition& pos);
+    int getMaterialID(const std::size_t element_id) const;
 
     Eigen::MatrixXd const& getPermeability(
+        const int material_id,
         const double t,
         const ProcessLib::SpatialPosition& pos,
         const int dim) const;
 
-    double getPorosity(const double t, const ProcessLib::SpatialPosition& pos,
-                       const double p, const double T,
-                       const double porosity_variable) const;
-
+    double getPorosity(const int material_id, const double t,
+                       const ProcessLib::SpatialPosition& pos, const double p,
+                       const double T, const double porosity_variable) const;
+    double getSaturation(const int material_id, const double t,
+                         const ProcessLib::SpatialPosition& pos, const double p,
+                         const double T, const double pc) const;
+    double getCapillaryPressure(const int material_id, const double t,
+                                const ProcessLib::SpatialPosition& pos,
+                                const double p, const double T,
+                                const double saturation) const;
+    double getSaturationDerivative(const int material_id, const double t,
+                                   const ProcessLib::SpatialPosition& pos,
+                                   const double p, const double T,
+                                   const double saturation) const;
     double getLiquidDensity(const double p, const double T) const;
     double getGasDensity(const double p, const double T) const;
     double getGasViscosity(const double p, const double T) const;
     double getLiquidViscosity(const double p, const double T) const;
-    double getDerivGasDensity(double const p, double const T) const;
-
-protected:
-    /// A flag to indicate whether the reference member, _material_ids,
-    /// is not assigned.
-    const bool _has_material_ids;
-
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _viscosity;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density;
-    std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity;
+    double getGasDensityDerivative(double const p, double const T) const;
+    double getNonwetRelativePermeability(const double t,
+                                         const ProcessLib::SpatialPosition& pos,
+                                         const double p, const double T,
+                                         const double saturation) const;
+    double getWetRelativePermeability(const double t,
+                                      const ProcessLib::SpatialPosition& pos,
+                                      const double p, const double T,
+                                      const double saturation) const;
 
+private:
     /** Use two phase models for different material zones.
-     *  Material IDs must be given as mesh element properties.
-     */
-    MeshLib::PropertyVector<int> const& _material_ids;
+    *  Material IDs must be given as mesh element properties.
+    */
+    boost::optional<MeshLib::PropertyVector<int> const&> const _material_ids;
+
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty> const _liquid_density;
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty> const _liquid_viscosity;
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty> const _gas_density;
+    std::unique_ptr<MaterialLib::Fluid::FluidProperty> const _gas_viscosity;
 
-    int _current_material_id = 0;
-    std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+    std::vector<Eigen::MatrixXd> const _intrinsic_permeability_models;
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> const
         _porosity_models;
-    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+    std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> const
         _storage_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>> const
+        _capillary_pressure_models;
+    std::vector<
+        std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>> const
+        _relative_permeability_models;
 };
 
 }  // end of namespace
-- 
GitLab