diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index cc74d66603d7530c12fb96643fb48cc5bebe65a9..ac841d3c3029c2d13f085efd90fe04262be4cc43 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -29,6 +29,41 @@ namespace ProcessLib
{
namespace HT
{
+
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> intrinsicPermeability(
+ MaterialPropertyLib::PropertyDataType const& values)
+{
+ if (boost::get<double>(&values))
+ {
+ return Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity() *
+ boost::get<double>(values);
+ }
+ if (boost::get<MaterialPropertyLib::Vector>(&values))
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, 1> const>(
+ boost::get<MaterialPropertyLib::Vector>(values).data(),
+ GlobalDim, 1)
+ .asDiagonal();
+ }
+ if (boost::get<MaterialPropertyLib::Tensor2d>(&values))
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, GlobalDim> const>(
+ boost::get<MaterialPropertyLib::Tensor2d>(values).data(), GlobalDim,
+ GlobalDim);
+ }
+ if (boost::get<MaterialPropertyLib::Tensor>(&values))
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, GlobalDim> const>(
+ boost::get<MaterialPropertyLib::Tensor>(values).data(), GlobalDim,
+ GlobalDim);
+ }
+ OGS_FATAL(
+ "Intrinsic permeability parameter values size is neither one nor %d "
+ "nor %d squared.",
+ GlobalDim, GlobalDim);
+}
+
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
class HTFEM : public HTLocalAssemblerInterface
@@ -93,7 +128,7 @@ public:
/// Computes the flux in the point \c pnt_local_coords that is given in
/// local coordinates using the values from \c local_x.
Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
- double const t,
+ double const /*t*/,
std::vector<double> const& local_x) const override
{
// eval dNdx and invJ at given point
@@ -128,16 +163,16 @@ public:
vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
p_int_pt;
- auto const K =
- this->_material_properties.porous_media_properties
- .getIntrinsicPermeability(t, pos)
- .getValue(t, pos, 0.0,
- boost::get<double>(vars[static_cast<int>(
- MaterialPropertyLib::Variable::temperature)]));
-
auto const& medium =
*_material_properties.media_map->getMedium(_element.getID());
auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
+
+ auto const K = intrinsicPermeability<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
+
auto const mu =
liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
.template value<double>(vars);
@@ -240,7 +275,7 @@ protected:
}
std::vector<double> const& getIntPtDarcyVelocityLocal(
- const double t, std::vector<double> const& local_p,
+ const double /*t*/, std::vector<double> const& local_p,
std::vector<double> const& local_T, std::vector<double>& cache) const
{
auto const n_integration_points =
@@ -262,6 +297,7 @@ protected:
auto const& medium =
*_material_properties.media_map->getMedium(_element.getID());
auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
@@ -281,9 +317,10 @@ protected:
vars[static_cast<int>(
MaterialPropertyLib::Variable::phase_pressure)] = p_int_pt;
- auto const K = _material_properties.porous_media_properties
- .getIntrinsicPermeability(t, pos)
- .getValue(t, pos, 0.0, T_int_pt);
+ auto const K = intrinsicPermeability<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
auto const mu =
liquid_phase
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 2c247adbd9fae49be27b5ac3fa91818a12ec6702..03ec0af2727fad7fd8eca2271329341a96ee6df3 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -101,6 +101,7 @@ public:
auto const& medium =
*process_data.media_map->getMedium(this->_element.getID());
auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
auto const& b = process_data.specific_body_force;
@@ -137,9 +138,13 @@ public:
auto const porosity =
process_data.porous_media_properties.getPorosity(t, pos)
.getValue(t, pos, 0.0, T_int_pt);
+
auto const intrinsic_permeability =
- process_data.porous_media_properties.getIntrinsicPermeability(
- t, pos).getValue(t, pos, 0.0, T_int_pt);
+ intrinsicPermeability<GlobalDim>(
+ solid_phase
+ .property(
+ MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
T_int_pt;
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index a625e0cb67d0d4ad029db2cbcda413b3adea4dc1..993e878f30956e61e2b6b46516b67561fd8c7184 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -128,10 +128,10 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
.getSpecificStorage(t, pos)
.getValue(0.0);
- auto const intrinsic_permeability =
- material_properties.porous_media_properties
- .getIntrinsicPermeability(t, pos)
- .getValue(t, pos, 0.0, T1_int_pt);
+ auto const intrinsic_permeability = intrinsicPermeability<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
// matrix assembly
@@ -208,6 +208,7 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
auto const& medium =
*material_properties.media_map->getMedium(this->_element.getID());
auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
auto const& b = material_properties.specific_body_force;
@@ -263,10 +264,10 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
.template value<double>(vars);
- auto const intrinsic_permeability =
- material_properties.porous_media_properties
- .getIntrinsicPermeability(t, pos)
- .getValue(t, pos, 0.0, T1_at_xi);
+ auto const intrinsic_permeability = intrinsicPermeability<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
GlobalDimVectorType const velocity =