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

RichardsMechanics refactoring step 3

See merge request ogs/ogs!5017

ProcessLib/RichardsMechanics/ConstitutiveRelations/ConstitutiveData.h

@@ -11,14 +11,19 @@
 #pragma once
 
 #include "Density.h"
+#include "DrySolidDensity.h"
 #include "LiquidDensity.h"
+#include "MicroPressure.h"
+#include "MicroSaturation.h"
 #include "ProcessLib/ConstitutiveRelations/Base.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Biot.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Saturation.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibilityData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h"
 #include "SaturationSecantDerivative.h"
 #include "StiffnessTensor.h"
@@ -30,15 +35,25 @@ template <int DisplacementDim>
 using StatefulData = std::tuple<
     StrainData<DisplacementDim>,
     ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
-        EffectiveStressData<DisplacementDim>>;
+        EffectiveStressData<DisplacementDim>,
+    ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
+        SwellingDataStateful<DisplacementDim>,
+    ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::
+        MechanicalStrainData<DisplacementDim>,
+    ProcessLib::ThermoRichardsMechanics::SaturationData,
+    ProcessLib::ThermoRichardsMechanics::PorosityData,
+    ProcessLib::ThermoRichardsMechanics::TransportPorosityData, MicroPressure,
+    MicroSaturation>;
 
 template <int DisplacementDim>
 using StatefulDataPrev = ProcessLib::ConstitutiveRelations::PrevStateOf<
     StatefulData<DisplacementDim>>;
 
 /// Data that is needed for output purposes, but not directly for the assembly.
 template <int DisplacementDim>
-using OutputData = std::tuple<>;
+using OutputData = std::tuple<
+    ProcessLib::ThermoRichardsMechanics::DarcyLawData<DisplacementDim>,
+    DrySolidDensity>;
 
 /// Data that is needed for the equation system assembly.
 template <int DisplacementDim>

ProcessLib/RichardsMechanics/ConstitutiveRelations/DrySolidDensity.h

@@ -0,0 +1,26 @@
+/**
+ * \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 <string_view>
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+// Apparent dry solid density
+using DrySolidDensity = BaseLib::StrongType<double, struct DrySolidDensityTag>;
+
+constexpr std::string_view ioName(struct DrySolidDensityTag*)
+{
+    return "dry_density_solid";
+}
+}  // namespace ProcessLib::RichardsMechanics

ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroPressure.h

@@ -0,0 +1,25 @@
+/**
+ * \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 <string_view>
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+using MicroPressure = BaseLib::StrongType<double, struct MicroPressureTag>;
+
+constexpr std::string_view ioName(struct MicroPressureTag*)
+{
+    return "micro_pressure";
+}
+}  // namespace ProcessLib::RichardsMechanics

ProcessLib/RichardsMechanics/ConstitutiveRelations/MicroSaturation.h

@@ -0,0 +1,25 @@
+/**
+ * \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 <string_view>
+
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::RichardsMechanics
+{
+using MicroSaturation = BaseLib::StrongType<double, struct MicroSaturationTag>;
+
+constexpr std::string_view ioName(struct MicroSaturationTag*)
+{
+    return "micro_saturation";
+}
+}  // namespace ProcessLib::RichardsMechanics

ProcessLib/RichardsMechanics/IntegrationPointData.h

@@ -12,7 +12,10 @@
 
 #include <memory>
 
+#include "ConstitutiveRelations/Base.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "MathLib/KelvinVector.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/SolidMechanics.h"
 
 namespace ProcessLib
 {
@@ -22,74 +25,24 @@ template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
           typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
 struct IntegrationPointData final
 {
-    explicit IntegrationPointData(
-        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
-            solid_material)
-        : solid_material(solid_material),
-          material_state_variables(
-              solid_material.createMaterialStateVariables())
-    {
-        // Initialize current time step values
-        static const int kelvin_vector_size =
-            MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-        sigma_sw.setZero(kelvin_vector_size);
-        eps_m.setZero(kelvin_vector_size);
-
-        // Previous time step values are not initialized and are set later.
-        eps_m_prev.resize(kelvin_vector_size);
-    }
-
-    typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
-    typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
-
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
     typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
 
     typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
     typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
 
-    typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy;
-
-    double liquid_pressure_m = std::numeric_limits<double>::quiet_NaN();
-    double liquid_pressure_m_prev = std::numeric_limits<double>::quiet_NaN();
-    double saturation = std::numeric_limits<double>::quiet_NaN();
-    double saturation_prev = std::numeric_limits<double>::quiet_NaN();
-    double saturation_m = std::numeric_limits<double>::quiet_NaN();
-    double saturation_m_prev = std::numeric_limits<double>::quiet_NaN();
-    double porosity = std::numeric_limits<double>::quiet_NaN();
-    double porosity_prev = std::numeric_limits<double>::quiet_NaN();
-    double transport_porosity = std::numeric_limits<double>::quiet_NaN();
-    double transport_porosity_prev = std::numeric_limits<double>::quiet_NaN();
-    double dry_density_solid = std::numeric_limits<double>::quiet_NaN();
-    double dry_density_pellet_saturated =
-        std::numeric_limits<double>::quiet_NaN();
-    double dry_density_pellet_unsaturated =
-        std::numeric_limits<double>::quiet_NaN();
-
-    MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
-    std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
-        DisplacementDim>::MaterialStateVariables>
-        material_state_variables;
     double integration_weight = std::numeric_limits<double>::quiet_NaN();
 
-    void pushBackState()
-    {
-        eps_m_prev = eps_m;
-        sigma_sw_prev = sigma_sw;
-        saturation_prev = saturation;
-        saturation_m_prev = saturation_m;
-        porosity_prev = porosity;
-        transport_porosity_prev = transport_porosity;
-        liquid_pressure_m_prev = liquid_pressure_m;
-        material_state_variables->pushBackState();
-    }
-
-    MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
-    computeElasticTangentStiffness(
-        double const t,
-        ParameterLib::SpatialPosition const& x_position,
-        double const dt,
-        double const temperature)
+    MathLib::KelvinVector::
+        KelvinMatrixType<DisplacementDim> static computeElasticTangentStiffness(
+            double const t,
+            ParameterLib::SpatialPosition const& x_position,
+            double const dt,
+            double const temperature,
+            MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+                solid_material,
+            typename MaterialLib::Solids::MechanicsBase<DisplacementDim>::
+                MaterialStateVariables const& material_state_variables)
     {
         namespace MPL = MaterialPropertyLib;
 
@@ -108,7 +61,7 @@ struct IntegrationPointData final
 
         auto&& solution = solid_material.integrateStress(
             variable_array_prev, variable_array, t, x_position, dt,
-            *material_state_variables);
+            material_state_variables);
 
         if (!solution)
         {
@@ -121,7 +74,7 @@ struct IntegrationPointData final
         return C;
     }
 
-    typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
+    static typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
         MaterialPropertyLib::VariableArray const& variable_array,
         double const t,
         ParameterLib::SpatialPosition const& x_position,
@@ -132,13 +85,25 @@ struct IntegrationPointData final
                 DisplacementDim>& sigma_eff,
         PrevState<ProcessLib::ThermoRichardsMechanics::
                       ConstitutiveStress_StrainTemperature::EffectiveStressData<
-                          DisplacementDim>> const& sigma_eff_prev)
+                          DisplacementDim>> const& sigma_eff_prev,
+        ProcessLib::ThermoRichardsMechanics::
+            ConstitutiveStress_StrainTemperature::MechanicalStrainData<
+                DisplacementDim> const&
+        /*eps_m*/,
+        PrevState<
+            ProcessLib::ThermoRichardsMechanics::
+                ConstitutiveStress_StrainTemperature::MechanicalStrainData<
+                    DisplacementDim>> const& eps_m_prev,
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            solid_material,
+        std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+            DisplacementDim>::MaterialStateVariables>& material_state_variables)
     {
         MaterialPropertyLib::VariableArray variable_array_prev;
         variable_array_prev.stress = sigma_eff_prev->sigma_eff;
         variable_array_prev.mechanical_strain
             .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
-                eps_m_prev);
+                eps_m_prev->eps_m);
         variable_array_prev.temperature = temperature;
 
         auto&& solution = solid_material.integrateStress(