diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index 92a8bd4cfa244568506e7ff36586f4a658a1326a..6af23f988eb9a916ddd63dc5d84d751fdb973d82 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -17,6 +17,7 @@ append_source_files(SOURCES Fluid/WaterVaporProperties)
append_source_files(SOURCES MPL)
append_source_files(SOURCES MPL/Properties)
append_source_files(SOURCES MPL/Components)
+append_source_files(SOURCES MPL/Utils)
append_source_files(SOURCES PorousMedium)
append_source_files(SOURCES PorousMedium/Porosity)
diff --git a/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.cpp b/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b0f1bb55aab13fdf2c11f06453c99c3baf24c42
--- /dev/null
+++ b/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.cpp
@@ -0,0 +1,39 @@
+/*
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on July 31, 2019, 12:10 PM
+ */
+
+#include "FormEffectiveThermalConductivity.h"
+
+#include "FormEigenTensor.h"
+
+namespace MaterialPropertyLib
+{
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> formEffectiveThermalConductivity(
+ MaterialPropertyLib::PropertyDataType const& solid_thermal_conductivity,
+ const double fluid_thermal_conductivity, const double porosity)
+{
+ return (1.0 - porosity) *
+ formEigenTensor<GlobalDim>(solid_thermal_conductivity) +
+ porosity * fluid_thermal_conductivity *
+ Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity();
+}
+
+template Eigen::Matrix<double, 1, 1> formEffectiveThermalConductivity<1>(
+ MaterialPropertyLib::PropertyDataType const& solid_thermal_conductivity,
+ const double fluid_thermal_conductivity, const double porosity);
+template Eigen::Matrix<double, 2, 2> formEffectiveThermalConductivity<2>(
+ MaterialPropertyLib::PropertyDataType const& solid_thermal_conductivity,
+ const double fluid_thermal_conductivity, const double porosity);
+template Eigen::Matrix<double, 3, 3> formEffectiveThermalConductivity<3>(
+ MaterialPropertyLib::PropertyDataType const& solid_thermal_conductivity,
+ const double fluid_thermal_conductivity, const double porosity);
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h b/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h
new file mode 100644
index 0000000000000000000000000000000000000000..52942f6ec65c219d6ee9f2e357d4909347bd7751
--- /dev/null
+++ b/MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h
@@ -0,0 +1,24 @@
+/*
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on July 31, 2019, 12:10 PM
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> formEffectiveThermalConductivity(
+ MaterialPropertyLib::PropertyDataType const& solid_thermal_conductivity,
+ const double fluid_thermal_conductivity, const double porosity);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Utils/FormEigenTensor.cpp b/MaterialLib/MPL/Utils/FormEigenTensor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6479e24fdbfa4b0898cd2936c39cc0354ad8343d
--- /dev/null
+++ b/MaterialLib/MPL/Utils/FormEigenTensor.cpp
@@ -0,0 +1,89 @@
+/*
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on July 31, 2019, 11:28 AM
+ */
+
+#include "FormEigenTensor.h"
+
+#include <boost/variant/static_visitor.hpp>
+
+#include "MaterialLib/MPL/PropertyType.h"
+
+namespace MaterialPropertyLib
+{
+template <int GlobalDim>
+struct FormEigenTensor
+ : boost::static_visitor<Eigen::Matrix<double, GlobalDim, GlobalDim>>
+{
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ double const& value) const
+ {
+ return Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity() * value;
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ MaterialPropertyLib::Vector const& values) const
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, 1> const>(
+ values.data(), GlobalDim, 1)
+ .asDiagonal();
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ MaterialPropertyLib::Tensor2d const& values) const
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, GlobalDim> const>(
+ values.data(), GlobalDim, GlobalDim);
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ MaterialPropertyLib::Tensor const& values) const
+ {
+ return Eigen::Map<Eigen::Matrix<double, GlobalDim, GlobalDim> const>(
+ values.data(), GlobalDim, GlobalDim);
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ MaterialPropertyLib::SymmTensor const& /*values*/) const
+ {
+ OGS_FATAL(
+ "The value of MaterialPropertyLib::SymmTensor is inapplicable");
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ std::string const& /*values*/) const
+ {
+ OGS_FATAL("The value of std::string is inapplicable");
+ }
+
+ Eigen::Matrix<double, GlobalDim, GlobalDim> operator()(
+ MaterialPropertyLib::Pair const& /*values*/) const
+ {
+ OGS_FATAL("The size of tensor is neither one nor %d nor %d squared.",
+ GlobalDim, GlobalDim);
+ }
+};
+
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> formEigenTensor(
+ MaterialPropertyLib::PropertyDataType const& values)
+{
+ return boost::apply_visitor(FormEigenTensor<GlobalDim>(), values);
+}
+
+template Eigen::Matrix<double, 1, 1> formEigenTensor<1>(
+ MaterialPropertyLib::PropertyDataType const& values);
+
+template Eigen::Matrix<double, 2, 2> formEigenTensor<2>(
+ MaterialPropertyLib::PropertyDataType const& values);
+
+template Eigen::Matrix<double, 3, 3> formEigenTensor<3>(
+ MaterialPropertyLib::PropertyDataType const& values);
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Utils/FormEigenTensor.h b/MaterialLib/MPL/Utils/FormEigenTensor.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c72b96de0f475e4f185ed4a6cdecfed347b5b2e
--- /dev/null
+++ b/MaterialLib/MPL/Utils/FormEigenTensor.h
@@ -0,0 +1,23 @@
+/*
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on July 29, 2019, 2:48 PM
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> formEigenTensor(
+ MaterialPropertyLib::PropertyDataType const& values);
+} // namespace MaterialPropertyLib
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index b85c26806427fc6d633cb816c3209ecf8bf5cf8b..27984d5db843b7c461e15a0bdae9ac61ed79614d 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -15,7 +15,8 @@
#include "HTMaterialProperties.h"
#include "MaterialLib/MPL/Medium.h"
-#include "MaterialLib/MPL/PropertyType.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -29,41 +30,6 @@ 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
@@ -165,7 +131,7 @@ public:
auto const& liquid_phase = medium.phase("AqueousLiquid");
auto const& solid_phase = medium.phase("Solid");
- auto const K = intrinsicPermeability<GlobalDim>(
+ auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
solid_phase
.property(MaterialPropertyLib::PropertyType::permeability)
.value(vars));
@@ -231,9 +197,8 @@ protected:
}
GlobalDimMatrixType getThermalConductivityDispersivity(
- MaterialPropertyLib::VariableArray const& vars,
- const double porosity, const double fluid_density,
- const double specific_heat_capacity_fluid,
+ MaterialPropertyLib::VariableArray const& vars, const double porosity,
+ const double fluid_density, const double specific_heat_capacity_fluid,
const GlobalDimVectorType& velocity, const GlobalDimMatrixType& I)
{
auto const& medium =
@@ -327,7 +292,7 @@ protected:
vars[static_cast<int>(
MaterialPropertyLib::Variable::phase_pressure)] = p_int_pt;
- auto const K = intrinsicPermeability<GlobalDim>(
+ auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
solid_phase
.property(MaterialPropertyLib::PropertyType::permeability)
.value(vars));
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index e86f9dc7da88bba5266a4ba4bd6034d2fbeafbc9..9471fcdc260080ae61cdf7417d31fa02ca332f8a 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -16,6 +16,8 @@
#include "HTMaterialProperties.h"
#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
@@ -139,7 +141,7 @@ public:
.template value<double>(vars);
auto const intrinsic_permeability =
- intrinsicPermeability<GlobalDim>(
+ MaterialPropertyLib::formEigenTensor<GlobalDim>(
solid_phase
.property(
MaterialPropertyLib::PropertyType::permeability)
@@ -162,7 +164,8 @@ public:
.template value<double>(vars);
// Use the viscosity model to compute the viscosity
- auto const viscosity = liquid_phase
+ auto const viscosity =
+ liquid_phase
.property(MaterialPropertyLib::PropertyType::viscosity)
.template value<double>(vars);
GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
diff --git a/ProcessLib/HT/StaggeredHTFEM-impl.h b/ProcessLib/HT/StaggeredHTFEM-impl.h
index 4254412621a8fa1f6f5c209afa48e23ec1a46507..90b380727ed715a13d62b1938c17d1d99cc86edc 100644
--- a/ProcessLib/HT/StaggeredHTFEM-impl.h
+++ b/ProcessLib/HT/StaggeredHTFEM-impl.h
@@ -14,7 +14,7 @@
#include "StaggeredHTFEM.h"
#include "MaterialLib/MPL/Medium.h"
-#include "MaterialLib/MPL/PropertyType.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
namespace ProcessLib
@@ -126,11 +126,13 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
solid_phase.property(MaterialPropertyLib::PropertyType::storage)
.template value<double>(vars);
- auto const intrinsic_permeability = intrinsicPermeability<GlobalDim>(
- solid_phase
- .property(MaterialPropertyLib::PropertyType::permeability)
- .value(vars));
- GlobalDimMatrixType const K_over_mu = intrinsic_permeability / viscosity;
+ auto const intrinsic_permeability =
+ MaterialPropertyLib::formEigenTensor<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
+ GlobalDimMatrixType const K_over_mu =
+ intrinsic_permeability / viscosity;
// matrix assembly
local_M.noalias() +=
@@ -261,10 +263,11 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
.template value<double>(vars);
- auto const intrinsic_permeability = intrinsicPermeability<GlobalDim>(
- solid_phase
- .property(MaterialPropertyLib::PropertyType::permeability)
- .value(vars));
+ auto const intrinsic_permeability =
+ MaterialPropertyLib::formEigenTensor<GlobalDim>(
+ solid_phase
+ .property(MaterialPropertyLib::PropertyType::permeability)
+ .value(vars));
GlobalDimMatrixType const K_over_mu =
intrinsic_permeability / viscosity;