[PL/TH2M] Stress/strain/saturation IP output.

Minimum set required for restart.

ProcessLib/TH2M/LocalAssemblerInterface.h

@@ -20,6 +20,10 @@ namespace TH2M
 struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
                                  public NumLib::ExtrapolatableElement
 {
+    virtual std::size_t setIPDataInitialConditions(
+        std::string const& name, double const* values,
+        int const integration_order) = 0;
+
     virtual std::vector<double> getSigma() const = 0;
 
     virtual std::vector<double> const& getIntPtSigma(
@@ -28,6 +32,8 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         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(
         const double t,
         std::vector<GlobalVector*> const& x,
@@ -76,6 +82,8 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> getSaturation() const = 0;
+
     virtual std::vector<double> const& getIntPtSaturation(
         const double t,
         std::vector<GlobalVector*> const& x,

ProcessLib/TH2M/TH2MFEM-impl.h

@@ -18,6 +18,7 @@
 #include "MathLib/KelvinVector.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 
 namespace ProcessLib
 {
@@ -341,6 +342,51 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     return ip_constitutive_variables;
 }
 
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+std::size_t TH2MLocalAssembler<
+    ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+    DisplacementDim>::setIPDataInitialConditions(std::string const& name,
+                                                 double const* values,
+                                                 int const integration_order)
+{
+    if (integration_order !=
+        static_cast<int>(_integration_method.getIntegrationOrder()))
+    {
+        OGS_FATAL(
+            "Setting integration point initial conditions; The integration "
+            "order of the local assembler for element {:d} is different "
+            "from the integration order in the initial condition.",
+            _element.getID());
+    }
+
+    if (name == "sigma_ip")
+    {
+        if (_process_data.initial_stress != nullptr)
+        {
+            OGS_FATAL(
+                "Setting initial conditions for stress from integration "
+                "point data and from a parameter '{:s}' is not possible "
+                "simultaneously.",
+                _process_data.initial_stress->name);
+        }
+        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+            values, _ip_data, &IpData::sigma_eff);
+    }
+
+    if (name == "saturation_ip")
+    {
+        return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
+                                                         &IpData::s_L);
+    }
+    if (name == "epsilon_ip")
+    {
+        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+            values, _ip_data, &IpData::eps);
+    }
+    return 0;
+}
+
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           typename IntegrationMethod, int DisplacementDim>
 void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,

ProcessLib/TH2M/TH2MFEM.h

@@ -28,6 +28,7 @@
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+#include "ProcessLib/Utils/TransposeInPlace.h"
 #include "TH2MProcessData.h"
 
 namespace ProcessLib
@@ -81,6 +82,12 @@ public:
                        TH2MProcessData<DisplacementDim>& process_data);
 
 private:
+    /// \return the number of read integration points.
+    std::size_t setIPDataInitialConditions(
+        std::string const& name,
+        double const* values,
+        int const integration_order) override;
+
     void setInitialConditionsConcrete(std::vector<double> const& local_x,
                                       double const t,
                                       bool const /*use_monolithic_scheme*/,
@@ -186,8 +193,15 @@ private:
     // There should be only one.
     std::vector<double> getSigma() const override
     {
-        return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
-            _ip_data, &IpData::sigma_eff);
+        {
+            constexpr int kelvin_vector_size =
+                MathLib::KelvinVector::kelvin_vector_dimensions(
+                    DisplacementDim);
+
+            return transposeInPlace<kelvin_vector_size>(
+                [this](std::vector<double>& values)
+                { return getIntPtSigma(0, {}, {}, values); });
+        }
     }
 
     std::vector<double> const& getIntPtSigma(
@@ -200,6 +214,16 @@ private:
             _ip_data, &IpData::sigma_eff, cache);
     }
 
+    std::vector<double> getEpsilon() const override
+    {
+        constexpr int kelvin_vector_size =
+            MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+
+        return transposeInPlace<kelvin_vector_size>(
+            [this](std::vector<double>& values)
+            { return getIntPtEpsilon(0, {}, {}, values); });
+    }
+
     virtual std::vector<double> const& getIntPtEpsilon(
         const double /*t*/,
         std::vector<GlobalVector*> const& /*x*/,
@@ -264,6 +288,13 @@ private:
                                                          cache);
     }
 
+    std::vector<double> getSaturation() const override
+    {
+        std::vector<double> result;
+        getIntPtSaturation(0, {}, {}, result);
+        return result;
+    }
+
     virtual std::vector<double> const& getIntPtSaturation(
         const double /*t*/,
         std::vector<GlobalVector*> const& /*x*/,

ProcessLib/TH2M/TH2MProcess.cpp

@@ -45,6 +45,70 @@ TH2MProcess<DisplacementDim>::TH2MProcess(
 
     _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
         mesh, "HydraulicFlow", MeshLib::MeshItemType::Node, 1);
+
+    // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
+    // properties, s.t. there is no "overlapping" with cell/point data.
+    // See getOrCreateMeshProperty.
+    _integration_point_writer.emplace_back(
+        std::make_unique<IntegrationPointWriter>(
+            "sigma_ip",
+            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+            integration_order,
+            [this]()
+            {
+                // Result containing integration point data for each local
+                // assembler.
+                std::vector<std::vector<double>> result;
+                result.resize(_local_assemblers.size());
+
+                for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+                {
+                    auto const& local_asm = *_local_assemblers[i];
+                    result[i] = local_asm.getSigma();
+                }
+
+                return result;
+            }));
+
+    _integration_point_writer.emplace_back(
+        std::make_unique<IntegrationPointWriter>(
+            "saturation_ip", 1 /*n components*/, integration_order,
+            [this]()
+            {
+                // Result containing integration point data for each local
+                // assembler.
+                std::vector<std::vector<double>> result;
+                result.resize(_local_assemblers.size());
+
+                for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+                {
+                    auto const& local_asm = *_local_assemblers[i];
+                    result[i] = local_asm.getSaturation();
+                }
+
+                return result;
+            }));
+
+    _integration_point_writer.emplace_back(
+        std::make_unique<IntegrationPointWriter>(
+            "epsilon_ip",
+            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+            integration_order,
+            [this]()
+            {
+                // Result containing integration point data for each local
+                // assembler.
+                std::vector<std::vector<double>> result;
+                result.resize(_local_assemblers.size());
+
+                for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+                {
+                    auto const& local_asm = *_local_assemblers[i];
+                    result[i] = local_asm.getEpsilon();
+                }
+
+                return result;
+            }));
 }
 
 template <int DisplacementDim>
@@ -234,6 +298,57 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
             const_cast<MeshLib::Mesh&>(mesh), "temperature_interpolated",
             MeshLib::MeshItemType::Node, 1);
 
+    // Set initial conditions for integration point data.
+    for (auto const& ip_writer : _integration_point_writer)
+    {
+        // Find the mesh property with integration point writer's name.
+        auto const& name = ip_writer->name();
+        if (!mesh.getProperties().existsPropertyVector<double>(name))
+        {
+            continue;
+        }
+        auto const& _meshproperty =
+            *mesh.getProperties().template getPropertyVector<double>(name);
+
+        // The mesh property must be defined on integration points.
+        if (_meshproperty.getMeshItemType() !=
+            MeshLib::MeshItemType::IntegrationPoint)
+        {
+            continue;
+        }
+
+        auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
+
+        // Check the number of components.
+        if (ip_meta_data.n_components !=
+            _meshproperty.getNumberOfGlobalComponents())
+        {
+            OGS_FATAL(
+                "Different number of components in meta data ({:d}) than in "
+                "the integration point field data for '{:s}': {:d}.",
+                ip_meta_data.n_components, name,
+                _meshproperty.getNumberOfGlobalComponents());
+        }
+
+        // Now we have a properly named vtk's field data array and the
+        // corresponding meta data.
+        std::size_t position = 0;
+        for (auto& local_asm : _local_assemblers)
+        {
+            std::size_t const integration_points_read =
+                local_asm->setIPDataInitialConditions(
+                    name, &_meshproperty[position],
+                    ip_meta_data.integration_order);
+            if (integration_points_read == 0)
+            {
+                OGS_FATAL(
+                    "No integration points read in the integration point "
+                    "initial conditions set function.");
+            }
+            position += integration_points_read * ip_meta_data.n_components;
+        }
+    }
+
     // Initialize local assemblers after all variables have been set.
     GlobalExecutor::executeMemberOnDereferenced(
         &LocalAssemblerInterface::initialize, _local_assemblers,