Merge branch 'rm-constitutive-setting-step-2' into 'master'

Further RichardsMechanics refactoring

See merge request ogs/ogs!5012

ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h

@@ -19,15 +19,18 @@
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h"
 #include "SaturationSecantDerivative.h"
 #include "StiffnessTensor.h"
 
 namespace ProcessLib::RichardsMechanics
 {
-// TODO directly declare these type aliases in Traits.h
 /// Data whose state must be tracked by the TRM process.
 template <int DisplacementDim>
-using StatefulData = std::tuple<>;
+using StatefulData = std::tuple<
+    StrainData<DisplacementDim>,
+    ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
+        EffectiveStressData<DisplacementDim>>;
 
 template <int DisplacementDim>
 using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf<

ProcessLib/RichardsMechanics/IntegrationPointData.h

@@ -32,20 +32,14 @@ struct IntegrationPointData final
         // Initialize current time step values
         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);
 
         // Previous time step values are not initialized and are set later.
-        eps_prev.resize(kelvin_vector_size);
         eps_m_prev.resize(kelvin_vector_size);
-        sigma_eff_prev.resize(kelvin_vector_size);
     }
 
-    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;
 
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
@@ -80,9 +74,7 @@ struct IntegrationPointData final
 
     void pushBackState()
     {
-        eps_prev = eps;
         eps_m_prev = eps_m;
-        sigma_eff_prev = sigma_eff;
         sigma_sw_prev = sigma_sw;
         saturation_prev = saturation;
         saturation_m_prev = saturation_m;
@@ -134,12 +126,16 @@ struct IntegrationPointData final
         double const t,
         ParameterLib::SpatialPosition const& x_position,
         double const dt,
-        double const temperature)
+        double const temperature,
+        ProcessLib::ThermoRichardsMechanics::
+            ConstitutiveStress_StrainTemperature::EffectiveStressData<
+                DisplacementDim>& sigma_eff,
+        PrevState<ProcessLib::ThermoRichardsMechanics::
+                      ConstitutiveStress_StrainTemperature::EffectiveStressData<
+                          DisplacementDim>> const& sigma_eff_prev)
     {
         MaterialPropertyLib::VariableArray variable_array_prev;
-        variable_array_prev.stress
-            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
-                sigma_eff_prev);
+        variable_array_prev.stress = sigma_eff_prev->sigma_eff;
         variable_array_prev.mechanical_strain
             .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                 eps_m_prev);
@@ -155,7 +151,8 @@ struct IntegrationPointData final
         }
 
         MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
-        std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+        std::tie(sigma_eff.sigma_eff, material_state_variables, C) =
+            std::move(*solution);
 
         return C;
     }

ProcessLib/RichardsMechanics/LocalAssemblerInterface.h

@@ -10,9 +10,14 @@
 
 #pragma once
 
+#include "ConstitutiveRelations/ConstitutiveData.h"
 #include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/MaterialState.h"
+#include "RichardsMechanicsProcessData.h"
 
 namespace ProcessLib
 {
@@ -22,6 +27,41 @@ template <int DisplacementDim>
 struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
                                  public NumLib::ExtrapolatableElement
 {
+    LocalAssemblerInterface(
+        MeshLib::Element const& e,
+        NumLib::GenericIntegrationMethod const& integration_method,
+        bool const is_axially_symmetric,
+        RichardsMechanicsProcessData<DisplacementDim>& process_data)
+        : process_data_(process_data),
+          integration_method_(integration_method),
+          element_(e),
+          is_axially_symmetric_(is_axially_symmetric)
+    {
+        unsigned const n_integration_points =
+            integration_method_.getNumberOfPoints();
+
+        current_states_.resize(n_integration_points);
+        prev_states_.resize(n_integration_points);
+        output_data_.resize(n_integration_points);
+
+        material_states_.reserve(n_integration_points);
+
+        auto const& solid_material =
+            MaterialLib::Solids::selectSolidConstitutiveRelation(
+                process_data_.solid_materials, process_data_.material_ids,
+                e.getID());
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            material_states_.emplace_back(
+                solid_material.createMaterialStateVariables());
+
+            // Set initial strain field to zero.
+            std::get<StrainData<DisplacementDim>>(current_states_[ip]).eps =
+                KelvinVector<DisplacementDim>::Zero();
+        }
+    }
+
     virtual std::size_t setIPDataInitialConditions(
         std::string_view const name, double const* values,
         int const integration_order) = 0;
@@ -116,7 +156,19 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
     virtual typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables const&
     getMaterialStateVariablesAt(unsigned /*integration_point*/) const = 0;
-};
 
+protected:
+    RichardsMechanicsProcessData<DisplacementDim>& process_data_;
+    NumLib::GenericIntegrationMethod const& integration_method_;
+    MeshLib::Element const& element_;
+    bool const is_axially_symmetric_;
+
+    std::vector<StatefulData<DisplacementDim>> current_states_;
+    std::vector<StatefulDataPrev<DisplacementDim>> prev_states_;
+    std::vector<
+        ProcessLib::ThermoRichardsMechanics::MaterialStateData<DisplacementDim>>
+        material_states_;
+    std::vector<OutputData<DisplacementDim>> output_data_;
+};
 }  // namespace RichardsMechanics
 }  // namespace ProcessLib

ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h

@@ -125,25 +125,29 @@ public:
     void initializeConcrete() override
     {
         unsigned const n_integration_points =
-            _integration_method.getNumberOfPoints();
+            this->integration_method_.getNumberOfPoints();
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
+            auto& SD = this->current_states_[ip];
             auto& ip_data = _ip_data[ip];
 
             ParameterLib::SpatialPosition const x_position{
-                std::nullopt, _element.getID(), ip,
-                MathLib::Point3d(
-                    NumLib::interpolateCoordinates<
-                        ShapeFunctionDisplacement,
-                        ShapeMatricesTypeDisplacement>(_element, ip_data.N_u))};
+                std::nullopt, this->element_.getID(), ip,
+                MathLib::Point3d(NumLib::interpolateCoordinates<
+                                 ShapeFunctionDisplacement,
+                                 ShapeMatricesTypeDisplacement>(this->element_,
+                                                                ip_data.N_u))};
 
             /// Set initial stress from parameter.
-            if (_process_data.initial_stress != nullptr)
+            if (this->process_data_.initial_stress != nullptr)
             {
-                ip_data.sigma_eff =
+                std::get<ProcessLib::ThermoRichardsMechanics::
+                             ConstitutiveStress_StrainTemperature::
+                                 EffectiveStressData<DisplacementDim>>(SD)
+                    .sigma_eff =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
-                        DisplacementDim>((*_process_data.initial_stress)(
+                        DisplacementDim>((*this->process_data_.initial_stress)(
                         std::numeric_limits<
                             double>::quiet_NaN() /* time independent */,
                         x_position));
@@ -154,6 +158,8 @@ public:
                 t, x_position, *ip_data.material_state_variables);
 
             ip_data.pushBackState();
+
+            this->prev_states_[ip] = SD;
         }
     }
 
@@ -163,12 +169,18 @@ public:
                               int const /*process_id*/) override
     {
         unsigned const n_integration_points =
-            _integration_method.getNumberOfPoints();
+            this->integration_method_.getNumberOfPoints();
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             _ip_data[ip].pushBackState();
         }
+
+        // TODO move to the local assembler interface
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            this->prev_states_[ip] = this->current_states_[ip];
+        }
     }
 
     void computeSecondaryVariableConcrete(
@@ -336,15 +348,12 @@ private:
         MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
         MPL::Medium const* const medium, TemperatureData const T_data,
         CapillaryPressureData<DisplacementDim> const& p_cap_data,
-        ConstitutiveData<DisplacementDim>& CD);
-
-    RichardsMechanicsProcessData<DisplacementDim>& _process_data;
+        ConstitutiveData<DisplacementDim>& CD,
+        StatefulData<DisplacementDim>& SD,
+        StatefulDataPrev<DisplacementDim> const& SD_prev);
 
     std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
 
-    NumLib::GenericIntegrationMethod const& _integration_method;
-    MeshLib::Element const& _element;
-    bool const _is_axially_symmetric;
     SecondaryData<
         typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
         _secondary_data;