[PL/TH2M] Collect thermal conductivities data

MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h

@@ -39,9 +39,9 @@ public:
                            double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,
                             Variable const variable,
-                            ParameterLib::SpatialPosition const& /*pos*/,
-                            double const /*t*/,
-                            double const /*dt*/) const override;
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t,
+                            double const dt) const override;
 
 private:
     ParameterLib::CoordinateSystem const* const local_coordinate_system_;

ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h

@@ -32,6 +32,7 @@
 #include "SolidMechanics.h"
 #include "SolidThermalExpansion.h"
 #include "Swelling.h"
+#include "ThermalConductivity.h"
 #include "TotalStress.h"
 #include "VapourPartialPressure.h"
 #include "Viscosity.h"
@@ -140,14 +141,12 @@ struct ConstitutiveTempData
     PorosityDerivativeData porosity_d_data;
     SolidDensityDerivativeData solid_density_d_data;
     SolidHeatCapacityData solid_heat_capacity_data;
+    ThermalConductivityData<DisplacementDim> thermal_conductivity_data;
 
     using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
     using DisplacementDimMatrix =
         Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
 
-    DisplacementDimMatrix dlambda_dp_GR;
-    DisplacementDimMatrix dlambda_dp_cap;
-    DisplacementDimMatrix dlambda_dT;
     DisplacementDimVector drho_GR_h_w_eff_dp_GR_Npart;
     DisplacementDimMatrix drho_GR_h_w_eff_dp_GR_gradNpart;
     DisplacementDimVector drho_LR_h_w_eff_dp_cap_Npart;

ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include "Base.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+template <int DisplacementDim>
+struct ThermalConductivityData
+{
+    GlobalDimMatrix<DisplacementDim> lambda;
+    // Currently unused, but there is a comment in TH2MFEM-impl.h referring to
+    // this matrix
+    // GlobalDimMatrix<DisplacementDim> dlambda_dp_GR;
+    GlobalDimMatrix<DisplacementDim> dlambda_dp_cap;
+    GlobalDimMatrix<DisplacementDim> dlambda_dT;
+};
+
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/IntegrationPointData.h

@@ -48,7 +48,6 @@ struct IntegrationPointData final
     double rho_u_eff = std::numeric_limits<double>::quiet_NaN();
     double rho_u_eff_prev = std::numeric_limits<double>::quiet_NaN();
 
-    GlobalDimMatrixType lambda;
     GlobalDimVectorType d_CG;
     GlobalDimVectorType d_WG;
     GlobalDimVectorType d_CL;

ProcessLib/TH2M/TH2MFEM-impl.h

@@ -273,11 +273,12 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         double const phi_S = 1. - ip_out.porosity_data.phi;
 
         // thermal conductivity
-        ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
-            medium
-                .property(
-                    MaterialPropertyLib::PropertyType::thermal_conductivity)
-                .value(vars, pos, t, dt));
+        ip_cv.thermal_conductivity_data.lambda =
+            MaterialPropertyLib::formEigenTensor<DisplacementDim>(
+                medium
+                    .property(
+                        MaterialPropertyLib::PropertyType::thermal_conductivity)
+                    .value(vars, pos, t, dt));
 
         ip_data.h_S = ip_cv.solid_heat_capacity_data() * T;
         auto const u_S = ip_data.h_S;
@@ -409,12 +410,12 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                           .dValue(vars, MPL::Variable::temperature, pos, t, dt))
                 : MPL::formEigenTensor<DisplacementDim>(0.);
 
-        ip_cv.dlambda_dp_cap =
+        ip_cv.thermal_conductivity_data.dlambda_dp_cap =
             dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
 
-        ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
-                           phi_S * dlambda_SR_dT -
-                           ip_cv.porosity_d_data.dphi_dT * lambdaSR;
+        ip_cv.thermal_conductivity_data.dlambda_dT =
+            phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
+            phi_S * dlambda_SR_dT - ip_cv.porosity_d_data.dphi_dT * lambdaSR;
 
         // From p_LR = p_GR - p_cap it follows for
         // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
@@ -1392,7 +1393,8 @@ void TH2MLocalAssembler<
 
         MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
 
-        KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
+        KTT.noalias() +=
+            gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
 
         fT.noalias() -= NTT * rho_u_eff_dot * w;
 
@@ -2053,7 +2055,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
                 temperature_index, temperature_index)
             .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
 
-        KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
+        KTT.noalias() +=
+            gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
 
         // d KTT/dp_GR * T
         // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
@@ -2073,13 +2076,16 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         local_Jac
             .template block<temperature_size, W_size>(temperature_index,
                                                       W_index)
-            .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
+            .noalias() += gradNTT *
+                          ip_cv.thermal_conductivity_data.dlambda_dp_cap *
+                          gradT * Np * w;
 
         // d KTT/dT * T
         local_Jac
             .template block<temperature_size, temperature_size>(
                 temperature_index, temperature_index)
-            .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
+            .noalias() += gradNTT * ip_cv.thermal_conductivity_data.dlambda_dT *
+                          gradT * NT * w;
 
         // fT_1
         fT.noalias() -= NTT * rho_u_eff_dot * w;