From 4d55379dc40a5de085dabdcbb78c0ec4ff148d18 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Wed, 10 Jan 2024 18:37:49 +0100
Subject: [PATCH] [PL/TH2M] Replace eps initialization, no ip data

Since eps became stateless, the initialization is also using
current displacement to compute the strains.
With this the epsilon_ip data is no longer necessary and does
not have to be reset for certain restart simulations.
---
 ProcessLib/TH2M/IntegrationPointData.h           |  3 ++-
 ProcessLib/TH2M/TH2MFEM-impl.h                   | 16 ++++++++++++++++
 .../square_1e1_2_matIDs_t_0.7000.vtu             |  1 -
 3 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index d96684e7400..a14e57ffeac 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -41,7 +41,8 @@ struct IntegrationPointData final
             MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
         sigma_eff.setZero(kelvin_vector_size);
         sigma_sw.setZero(kelvin_vector_size);
-        eps.setZero(kelvin_vector_size);
+        eps.resize(
+            kelvin_vector_size);  // Later initialization from displacement
         eps_m.setZero(kelvin_vector_size);
         eps_m_prev.resize(kelvin_vector_size);
         sigma_eff_prev.resize(kelvin_vector_size);
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index f0f2fd2c691..9dcabc5d0a1 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -947,6 +947,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     auto const temperature =
         local_x.template segment<temperature_size>(temperature_index);
 
+    auto const displacement =
+        local_x.template segment<displacement_size>(displacement_index);
+
     constexpr double dt = std::numeric_limits<double>::quiet_NaN();
     auto const& medium = *_process_data.media_map.getMedium(_element.getID());
     auto const& solid_phase = medium.phase("Solid");
@@ -961,6 +964,12 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         auto& ip_data = _ip_data[ip];
         auto const& Np = ip_data.N_p;
         auto const& NT = Np;
+        auto const& Nu = ip_data.N_u;
+        auto const& gradNu = ip_data.dNdx_u;
+        auto const x_coord =
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, Nu);
         ParameterLib::SpatialPosition const pos{
             std::nullopt, _element.getID(), ip,
             MathLib::Point3d(
@@ -974,7 +983,14 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         double const T = NT.dot(temperature);
         vars.temperature = T;
 
+        auto const Bu =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                gradNu, Nu, x_coord, _is_axially_symmetric);
+
         auto& eps = ip_data.eps;
+        eps.noalias() = Bu * displacement;
 
         // Set volumetric strain rate for the general case without swelling.
         vars.volumetric_strain = Invariants::trace(eps);
diff --git a/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu b/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu
index e189ca76043..b94e9cc3c03 100644
--- a/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu
+++ b/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu
@@ -4,7 +4,6 @@
     <FieldData>
       <DataArray type="Int8" Name="IntegrationPointMetaData" NumberOfTuples="1002" format="appended" RangeMin="34"                   RangeMax="125"                  offset="0"                   />
       <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="25" format="appended" RangeMin="45"                   RangeMax="121"                  offset="312"                 />
-      <DataArray type="Float64" Name="epsilon_ip" NumberOfComponents="4" NumberOfTuples="48" format="appended" RangeMin="0.00049497474683"     RangeMax="0.00049497474683"     offset="400"                 />
       <DataArray type="Float64" Name="material_state_variable_ElasticStrain_ip" NumberOfComponents="4" NumberOfTuples="24" format="appended" RangeMin="0.00043074710448"     RangeMax="0.00043074710448"     offset="1000"                />
       <DataArray type="Float64" Name="material_state_variable_EquivalentPlasticStrain_ip" NumberOfTuples="24" format="appended" RangeMin="9.19908e-05"          RangeMax="9.19908e-05"          offset="1576"                />
       <DataArray type="Float64" Name="material_state_variable_damage.kappa_d_ip" NumberOfTuples="24" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1728"                />