From c64f0a1d9000db812509516ee6c8d64b9c5c71da Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Wed, 6 Apr 2022 17:03:03 +0200
Subject: [PATCH] [PL/TH2M] Add swelling stress ip variable

Adding output and integration point writer of yet unused
and uninitialized swelling stress variable.
---
 ProcessLib/TH2M/IntegrationPointData.h    |  3 +++
 ProcessLib/TH2M/LocalAssemblerInterface.h |  8 ++++++++
 ProcessLib/TH2M/TH2MFEM-impl.h            |  5 +++++
 ProcessLib/TH2M/TH2MFEM.h                 | 23 +++++++++++++++++++++++
 ProcessLib/TH2M/TH2MProcess.cpp           | 11 +++++++++++
 5 files changed, 50 insertions(+)

diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 0ec28659c41..85c00b68924 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -40,6 +40,7 @@ struct IntegrationPointData final
         static const int kelvin_vector_size =
             MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
         sigma_eff.setZero(kelvin_vector_size);
+        sigma_sw.setZero(kelvin_vector_size);
         eps.setZero(kelvin_vector_size);
         eps_m.setZero(kelvin_vector_size);
         eps_m_prev.resize(kelvin_vector_size);
@@ -53,6 +54,7 @@ struct IntegrationPointData final
         DisplacementDim, NPoints * DisplacementDim>
         N_u_op;
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
+    typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
     typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
 
@@ -183,6 +185,7 @@ struct IntegrationPointData final
         eps_prev = eps;
         eps_m_prev = eps_m;
         sigma_eff_prev = sigma_eff;
+        sigma_sw_prev = sigma_sw;
         s_L_prev = s_L;
 
         rho_G_h_G_prev = rho_G_h_G;
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index b66e9f7540f..53a5b52df1d 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -32,6 +32,14 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> getSwellingStress() const = 0;
+
+    virtual std::vector<double> const& getIntPtSwellingStress(
+        const double t,
+        std::vector<GlobalVector*> const& x,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+        std::vector<double>& cache) const = 0;
+
     virtual std::vector<double> getEpsilon() const = 0;
 
     virtual std::vector<double> const& getIntPtEpsilon(
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index f2885f03a32..b68fc1f837b 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -815,6 +815,11 @@ std::size_t TH2MLocalAssembler<
         return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                          &IpData::s_L);
     }
+    if (name == "swelling_stress_ip")
+    {
+        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+            values, _ip_data, &IpData::sigma_sw);
+    }
     if (name == "epsilon_ip")
     {
         return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 65b758f8490..1cb48cb87c1 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -215,6 +215,29 @@ private:
             _ip_data, &IpData::sigma_eff, cache);
     }
 
+    std::vector<double> getSwellingStress() const override
+    {
+        {
+            constexpr int kelvin_vector_size =
+                MathLib::KelvinVector::kelvin_vector_dimensions(
+                    DisplacementDim);
+
+            return transposeInPlace<kelvin_vector_size>(
+                [this](std::vector<double>& values)
+                { return getIntPtSwellingStress(0, {}, {}, values); });
+        }
+    }
+
+    std::vector<double> const& getIntPtSwellingStress(
+        const double /*t*/,
+        std::vector<GlobalVector*> const& /*x*/,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+        std::vector<double>& cache) const override
+    {
+        return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+            _ip_data, &IpData::sigma_sw, cache);
+    }
+
     std::vector<double> getEpsilon() const override
     {
         constexpr int kelvin_vector_size =
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index 2084a2c522c..1ac1744d9f1 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -56,6 +56,13 @@ TH2MProcess<DisplacementDim>::TH2MProcess(
             integration_order, _local_assemblers,
             &LocalAssemblerInterface::getSigma));
 
+    _integration_point_writer.emplace_back(
+        std::make_unique<IntegrationPointWriter>(
+            "swelling_stress_ip",
+            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+            integration_order, _local_assemblers,
+            &LocalAssemblerInterface::getSwellingStress));
+
     _integration_point_writer.emplace_back(
         std::make_unique<IntegrationPointWriter>(
             "saturation_ip", 1 /*n components*/, integration_order,
@@ -177,6 +184,10 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
                            MathLib::KelvinVector::KelvinVectorType<
                                DisplacementDim>::RowsAtCompileTime,
                            &LocalAssemblerInterface::getIntPtSigma);
+    add_secondary_variable("swelling_stress",
+                           MathLib::KelvinVector::KelvinVectorType<
+                               DisplacementDim>::RowsAtCompileTime,
+                           &LocalAssemblerInterface::getIntPtSwellingStress);
     add_secondary_variable("epsilon",
                            MathLib::KelvinVector::KelvinVectorType<
                                DisplacementDim>::RowsAtCompileTime,
-- 
GitLab