From 56d4677f9b65fcb57028fce7ab08d5bcff8928a4 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Wed, 17 Apr 2024 19:08:58 +0200
Subject: [PATCH] [PL/TH2M] Extract KuT part of U1 equation

---
 .../ConstitutiveRelations/ConstitutiveData.h  |  1 +
 .../ConstitutiveModels.h                      |  1 +
 .../TH2M/ConstitutiveRelations/UEquation.cpp  | 13 ++++++++++++
 .../TH2M/ConstitutiveRelations/UEquation.h    | 20 +++++++++++++++++++
 ProcessLib/TH2M/TH2MFEM-impl.h                | 15 ++++++++------
 ProcessLib/TH2M/TH2MFEM.h                     |  3 +++
 6 files changed, 47 insertions(+), 6 deletions(-)

diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index fa7c57299a8..38f09a60085 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -216,6 +216,7 @@ struct DerivativesData
     FW4LWpCDerivativeData<DisplacementDim> dfW_4_LWpC;
     FT1DerivativeData dfT_1;
     FT2DerivativeData<DisplacementDim> dfT_2;
+    FU1KUTDerivativeData<DisplacementDim> dfu_1_KuT;
     FU2KUpCDerivativeData dfu_2_KupC;
 };
 
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index c79d9e81fa8..277093d023f 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -116,6 +116,7 @@ struct ConstitutiveModels
     FT2Model<DisplacementDim> fT_2_model;
     FT3Model<DisplacementDim> fT_3_model;
 
+    FU1KUTModel<DisplacementDim> fu_1_KuT_model;
     FU2KUpCModel fu_2_KupC_model;
 };
 }  // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
index 07f391f74d8..dcea844e0a4 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
@@ -13,6 +13,19 @@ namespace ProcessLib::TH2M
 {
 namespace ConstitutiveRelations
 {
+template <int DisplacementDim>
+void FU1KUTModel<DisplacementDim>::dEval(
+    SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+    SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+    FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const
+{
+    dfu_1_KuT.dT = s_mech_data.stiffness_tensor *
+                   s_therm_exp_data.solid_linear_thermal_expansivity_vector;
+}
+
+template struct FU1KUTModel<2>;
+template struct FU1KUTModel<3>;
+
 void FU2KUpCModel::eval(BiotData const& biot_data,
                         BishopsData const& chi_S_L,
                         FU2KUpCData& fu_2_KupC) const
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
index bffa79f508b..70388bf50d2 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
@@ -13,11 +13,31 @@
 #include "Biot.h"
 #include "Bishops.h"
 #include "Saturation.h"
+#include "SolidMechanics.h"
+#include "SolidThermalExpansion.h"
 
 namespace ProcessLib::TH2M
 {
 namespace ConstitutiveRelations
 {
+template <int DisplacementDim>
+struct FU1KUTDerivativeData
+{
+    KelvinVector<DisplacementDim> dT;
+};
+
+template <int DisplacementDim>
+struct FU1KUTModel
+{
+    void dEval(
+        SolidMechanicsDataStateless<DisplacementDim> const& s_mech_data,
+        SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+        FU1KUTDerivativeData<DisplacementDim>& dfu_1_KuT) const;
+};
+
+extern template struct FU1KUTModel<2>;
+extern template struct FU1KUTModel<3>;
+
 struct FU2KUpCData
 {
     double m = nan;
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 3e0beb96ac3..a571f7ff8ea 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -489,6 +489,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     updateConstitutiveVariablesDerivatives(
         Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
         double const t, double const dt,
+        std::vector<
+            ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+            ip_constitutive_data,
         std::vector<
             ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
             ip_constitutive_variables,
@@ -524,6 +527,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     {
         auto const& ip_data = _ip_data[ip];
         auto& ip_dd = ip_d_data[ip];
+        auto const& ip_cd = ip_constitutive_data[ip];
         auto const& ip_cv = ip_constitutive_variables[ip];
         auto const& ip_out = this->output_data_[ip];
         auto const& current_state = this->current_states_[ip];
@@ -714,6 +718,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             ip_cv.viscosity_data,
             ip_dd.dfT_2);
 
+        models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
+                                    ip_dd.dfu_1_KuT);
+
         models.fu_2_KupC_model.dEval(ip_cv.biot_data,
                                      ip_cv.chi_S_L,
                                      pCap_data,
@@ -1339,7 +1346,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
         Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
                                           local_x_prev.size()),
-        t, dt, ip_constitutive_variables, models);
+        t, dt, ip_constitutive_data, ip_constitutive_variables, models);
 
     for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
     {
@@ -1781,11 +1788,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         local_Jac
             .template block<displacement_size, temperature_size>(
                 displacement_index, temperature_index)
-            .noalias() -=
-            BuT *
-            (ip_cd.s_mech_data.stiffness_tensor *
-             ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) *
-            NT * w;
+            .noalias() -= BuT * ip_dd.dfu_1_KuT.dT * NT * w;
 
         /* TODO (naumov) Test with gravity needed to check this Jacobian part.
         local_Jac
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 1ea5fcb3559..bad9ef32fd5 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -202,6 +202,9 @@ private:
     updateConstitutiveVariablesDerivatives(
         Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
         double const t, double const dt,
+        std::vector<
+            ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
+            ip_constitutive_data,
         std::vector<
             ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
             ip_constitutive_variables,
-- 
GitLab