diff --git a/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp index 9baa3b86db520af7a18a2622bb2e5300500e451e..2d1a496e6b952d4fb7676b1667873322735eb746 100644 --- a/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp +++ b/ProcessLib/ThermoRichardsMechanics/CreateThermoRichardsMechanicsProcess.cpp @@ -187,6 +187,11 @@ std::unique_ptr<Process> createThermoRichardsMechanicsProcessStage2( config.getConfigParameter<bool>("apply_body_force_for_deformation", true); + bool const initialize_porosity_from_medium_property = + //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_MECHANICS__initialize_porosity_from_medium_property} + config.getConfigParameter("initialize_porosity_from_medium_property", + true); + const bool use_TaylorHood_elements = variable_p->getShapeFunctionOrder() != variable_u->getShapeFunctionOrder() @@ -201,7 +206,9 @@ std::unique_ptr<Process> createThermoRichardsMechanicsProcessStage2( specific_body_force, mass_lumping, use_TaylorHood_elements, - apply_body_force_for_deformation}; + apply_body_force_for_deformation, + InitializePorosityFromMediumProperty{ + initialize_porosity_from_medium_property}}; SecondaryVariableCollection secondary_variables; diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h index 25f220b15ce8617bd0c78a20e12b12a8c6ee9d6c..51fa592a83adcad61e77c880b0a3dcdb8bfe18e1 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM.h @@ -299,23 +299,27 @@ public: time_independent, x_position)); } - // Initial porosity. Could be read from integration point data - // or mesh. - current_state.poro_data.phi = medium.property(MPL::porosity) - .template initialValue<double>( - x_position, time_independent); - - if (medium.hasProperty(MPL::PropertyType::transport_porosity)) + if (*this->process_data_.initialize_porosity_from_medium_property) { - current_state.transport_poro_data.phi = - medium.property(MPL::transport_porosity) + // Initial porosity. Could be read from integration point data + // or mesh. + current_state.poro_data.phi = + medium.property(MPL::porosity) .template initialValue<double>(x_position, time_independent); - } - else - { - current_state.transport_poro_data.phi = - current_state.poro_data.phi; + + if (medium.hasProperty(MPL::PropertyType::transport_porosity)) + { + current_state.transport_poro_data.phi = + medium.property(MPL::transport_porosity) + .template initialValue<double>(x_position, + time_independent); + } + else + { + current_state.transport_poro_data.phi = + current_state.poro_data.phi; + } } double const t = 0; // TODO (naumov) pass t from top diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h index 3391531c1d81770e17cf12979c8e1db2edd144b0..3c3958749d1c83aaef7c9d5771f2ec3b03b1c398 100644 --- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h +++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcessData.h @@ -14,6 +14,7 @@ #include <memory> #include <utility> +#include "BaseLib/StrongType.h" #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h" #include "ParameterLib/Parameter.h" @@ -21,6 +22,10 @@ namespace ProcessLib { namespace ThermoRichardsMechanics { +/// A type helping to avoid confusion of different boolean values. +using InitializePorosityFromMediumProperty = + BaseLib::StrongType<bool, struct InitializePorosityFromMediumPropertyTag>; + template <int DisplacementDim, typename ConstitutiveTraits> struct ThermoRichardsMechanicsProcessData { @@ -48,6 +53,9 @@ struct ThermoRichardsMechanicsProcessData bool const apply_body_force_for_deformation; + InitializePorosityFromMediumProperty const + initialize_porosity_from_medium_property; + MeshLib::PropertyVector<double>* element_saturation = nullptr; MeshLib::PropertyVector<double>* element_porosity = nullptr; MeshLib::PropertyVector<double>* element_liquid_density = nullptr;