diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index e3bc71e2a975f1fd12af67fe674e5ac672662867..cc41098b64d79e3d0df01441ad92cb9fc25e6859 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -178,8 +178,10 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
pressure_size, displacement_size>::Zero(pressure_size,
displacement_size);
- MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material =
- *_process_data.solid_materials[0];
+ auto const& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials, _process_data.material_ids,
+ _element.getID());
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -536,8 +538,10 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
pressure_size);
- MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material =
- *_process_data.solid_materials[0];
+ auto const& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials, _process_data.material_ids,
+ _element.getID());
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
index 683b8e619c9f2aa77cc2b779ce1e5f21c8bceb4e..1f4e41b8d3ae6d9dd26b2d68d9835a444775e5fa 100644
--- a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
+++ b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
@@ -72,6 +72,51 @@ struct IntegrationPointData final
material_state_variables->pushBackState();
}
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
+ computeElasticTangentStiffness(
+ double const t,
+ ParameterLib::SpatialPosition const& x_position,
+ double const dt,
+ double const temperature)
+ {
+ namespace MPL = MaterialPropertyLib;
+
+ MPL::VariableArray variable_array;
+ MPL::VariableArray variable_array_prev;
+
+ auto const null_state = solid_material.createMaterialStateVariables();
+
+ using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+
+ variable_array[static_cast<int>(MPL::Variable::stress)].emplace<KV>(
+ KV::Zero());
+ variable_array[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<KV>(KV::Zero());
+ variable_array[static_cast<int>(MPL::Variable::temperature)]
+ .emplace<double>(temperature);
+
+ variable_array_prev[static_cast<int>(MPL::Variable::stress)]
+ .emplace<KV>(KV::Zero());
+ variable_array_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<KV>(KV::Zero());
+ variable_array_prev[static_cast<int>(MPL::Variable::temperature)]
+ .emplace<double>(temperature);
+
+ auto&& solution =
+ solid_material.integrateStress(variable_array_prev, variable_array,
+ t, x_position, dt, *null_state);
+
+ if (!solution)
+ {
+ OGS_FATAL("Computation of elastic tangent stiffness failed.");
+ }
+
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C =
+ std::move(std::get<2>(*solution));
+
+ return C;
+ }
+
typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
MaterialPropertyLib::VariableArray const& variable_array,
double const t,
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index a82fe3d8bd352bd8866158283e27889778fefdef..7ca70890bb58fee5b2505c37309ede63d211312e 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -224,8 +224,10 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
KuT;
KuT.setZero(displacement_size, temperature_size);
- MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material =
- *_process_data.solid_materials[0];
+ auto const& solid_material =
+ MaterialLib::Solids::selectSolidConstitutiveRelation(
+ _process_data.solid_materials, _process_data.material_ids,
+ _element.getID());
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -297,8 +299,10 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
->property(MaterialPropertyLib::PropertyType::biot_coefficient)
.template value<double>(vars, x_position, t, dt);
+ auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
+ t, x_position, dt, static_cast<double>(T_int_pt));
auto const solid_skeleton_compressibility =
- 1 / solid_material.getBulkModulus(t, x_position);
+ 1 / solid_material.getBulkModulus(t, x_position, &C_el);
auto const beta_SR = (1 - alpha) * solid_skeleton_compressibility;
diff --git a/Tests/Data/ThermoHydroMechanics/BGRaCreepAndInitialStressAtIP_AREHS/arehs-salt-THM0_ts_5_t_157680000000.000000.vtu b/Tests/Data/ThermoHydroMechanics/BGRaCreepAndInitialStressAtIP_AREHS/arehs-salt-THM0_ts_5_t_157680000000.000000.vtu
index 2154671f9d9225ac9b21a7cd9599d5d20ec5b7d8..a4d4c862412e66e1711a974c7d03a2ff44a019d9 100644
--- a/Tests/Data/ThermoHydroMechanics/BGRaCreepAndInitialStressAtIP_AREHS/arehs-salt-THM0_ts_5_t_157680000000.000000.vtu
+++ b/Tests/Data/ThermoHydroMechanics/BGRaCreepAndInitialStressAtIP_AREHS/arehs-salt-THM0_ts_5_t_157680000000.000000.vtu
@@ -3,41 +3,42 @@
<UnstructuredGrid>
<FieldData>
<DataArray type="Int8" Name="IntegrationPointMetaData" NumberOfTuples="166" format="appended" RangeMin="34" RangeMax="125" offset="0" />
- <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="22" format="appended" RangeMin="45" RangeMax="103" offset="176" />
- <DataArray type="Float64" Name="epsilon_ip" NumberOfComponents="4" NumberOfTuples="53144" format="appended" RangeMin="3.4763901153e-05" RangeMax="0.025479754859" offset="260" />
- <DataArray type="Float64" Name="sigma_ip" NumberOfComponents="4" NumberOfTuples="53144" format="appended" RangeMin="66461.401901" RangeMax="97443995.414" offset="1668532" />
+ <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="26" format="appended" RangeMin="45" RangeMax="121" offset="176" />
+ <DataArray type="Float64" Name="epsilon_ip" NumberOfComponents="4" NumberOfTuples="53144" format="appended" RangeMin="3.4763898599e-05" RangeMax="0.025479754977" offset="268" />
+ <DataArray type="Float64" Name="sigma_ip" NumberOfComponents="4" NumberOfTuples="53144" format="appended" RangeMin="66461.443003" RangeMax="97443992.737" offset="1671740" />
</FieldData>
<Piece NumberOfPoints="15413" NumberOfCells="7592" >
<PointData>
- <DataArray type="Float64" Name="HydraulicFlow" format="appended" RangeMin="-642.2503427" RangeMax="249.55970989" offset="3819476" />
- <DataArray type="Float64" Name="MaterialForces" NumberOfComponents="2" format="appended" RangeMin="0.23136438348" RangeMax="7044778.7867" offset="3846352" />
- <DataArray type="Float64" Name="NodalForces" NumberOfComponents="2" format="appended" RangeMin="3.1622776602e+149" RangeMax="-nan" offset="4146556" />
- <DataArray type="Float64" Name="displacement" NumberOfComponents="2" format="appended" RangeMin="0" RangeMax="2.5757284198" offset="4172804" />
- <DataArray type="Float64" Name="epsilon" NumberOfComponents="4" format="appended" RangeMin="3.4758153185e-05" RangeMax="0.017045405211" offset="4484916" />
- <DataArray type="Float64" Name="pressure" format="appended" RangeMin="-2396594.8722" RangeMax="22354705.149" offset="4980492" />
- <DataArray type="Float64" Name="pressure_interpolated" format="appended" RangeMin="-2396594.8722" RangeMax="22354705.149" offset="5020296" />
- <DataArray type="Float64" Name="sigma" NumberOfComponents="4" format="appended" RangeMin="72591.824801" RangeMax="87725771.336" offset="5176040" />
- <DataArray type="Float64" Name="temperature" format="appended" RangeMin="0" RangeMax="365.64628245" offset="5802488" />
- <DataArray type="Float64" Name="temperature_interpolated" format="appended" RangeMin="279.42790834" RangeMax="365.64628245" offset="5839704" />
+ <DataArray type="Float64" Name="HeatFlux" format="appended" RangeMin="-642.25033413" RangeMax="249.55970685" offset="3822668" />
+ <DataArray type="Float64" Name="HydraulicFlow" format="appended" RangeMin="-0.0004782915916" RangeMax="0.00046890330617" offset="3850864" />
+ <DataArray type="Float64" Name="MaterialForces" NumberOfComponents="2" format="appended" RangeMin="0.23136438348" RangeMax="7044778.7867" offset="3876684" />
+ <DataArray type="Float64" Name="NodalForces" NumberOfComponents="2" format="appended" RangeMin="5.9604644775e-08" RangeMax="3789150972.1" offset="4176888" />
+ <DataArray type="Float64" Name="displacement" NumberOfComponents="2" format="appended" RangeMin="0" RangeMax="2.57572842" offset="4365100" />
+ <DataArray type="Float64" Name="epsilon" NumberOfComponents="4" format="appended" RangeMin="3.4758150612e-05" RangeMax="0.017045405274" offset="4677240" />
+ <DataArray type="Float64" Name="pressure" format="appended" RangeMin="-2396594.8902" RangeMax="22354704.523" offset="5172864" />
+ <DataArray type="Float64" Name="pressure_interpolated" format="appended" RangeMin="-2396594.8902" RangeMax="22354704.523" offset="5212664" />
+ <DataArray type="Float64" Name="sigma" NumberOfComponents="4" format="appended" RangeMin="72591.856344" RangeMax="87725769.416" offset="5368396" />
+ <DataArray type="Float64" Name="temperature" format="appended" RangeMin="0" RangeMax="365.64628247" offset="5994760" />
+ <DataArray type="Float64" Name="temperature_interpolated" format="appended" RangeMin="279.42790834" RangeMax="365.64628247" offset="6031992" />
</PointData>
<CellData>
- <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0" RangeMax="3" offset="5988164" />
- <DataArray type="Float64" Name="principal_stress_values" NumberOfComponents="3" format="appended" RangeMin="378818.44382" RangeMax="86370314.336" offset="5988292" />
- <DataArray type="Float64" Name="principal_stress_vector_1" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6217148" />
- <DataArray type="Float64" Name="principal_stress_vector_2" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6363504" />
- <DataArray type="Float64" Name="principal_stress_vector_3" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6425232" />
+ <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0" RangeMax="3" offset="6180460" />
+ <DataArray type="Float64" Name="principal_stress_values" NumberOfComponents="3" format="appended" RangeMin="378818.44382" RangeMax="86370314.336" offset="6180588" />
+ <DataArray type="Float64" Name="principal_stress_vector_1" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6409444" />
+ <DataArray type="Float64" Name="principal_stress_vector_2" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6555800" />
+ <DataArray type="Float64" Name="principal_stress_vector_3" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1" offset="6617528" />
</CellData>
<Points>
- <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="6708.2039325" offset="6525764" />
+ <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="6708.2039325" offset="6718060" />
</Points>
<Cells>
- <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="6838512" />
- <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="6975980" />
- <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="6990752" />
+ <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="7030808" />
+ <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="7168276" />
+ <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="7183048" />
</Cells>
</Piece>
</UnstructuredGrid>
<AppendedData encoding="base64">