[PL] SD: Set integration point initial conditions.

ProcessLib/Output/IntegrationPointWriter.cpp

@@ -71,6 +71,22 @@ static void addIntegrationPointMetaData(
               std::back_inserter(dictionary));
 }
 
+/// For the given json object and the name extract integration point meta data,
+/// or fail if no meta data was found for the given name.
+static ProcessLib::IntegrationPointMetaData extractIntegrationPointMetaData(
+    json const& meta_data, std::string const& name)
+{
+    for (auto const& md : meta_data["integration_point_arrays"])
+    {
+        if (md["name"] == name)
+        {
+            return {name, md["number_of_components"], md["integration_order"]};
+        }
+    }
+    OGS_FATAL("No integration point meta data with name \"%s\" found.",
+              name.c_str());
+}
+
 namespace ProcessLib
 {
 void addIntegrationPointWriter(
@@ -88,4 +104,36 @@ void addIntegrationPointWriter(
         addIntegrationPointMetaData(mesh, meta_data);
     }
 }
+
+IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
+                                                     std::string const& name)
+{
+    if (!mesh.getProperties().existsPropertyVector<char>(
+            "IntegrationPointMetaData"))
+    {
+        OGS_FATAL(
+            "Integration point data '%s' is present in the vtk field "
+            "data but the required \"IntegrationPointMetaData\" array "
+            "is not available.",
+            name.c_str());
+    }
+    auto const& mesh_property_ip_meta_data =
+        *mesh.getProperties().template getPropertyVector<char>(
+            "IntegrationPointMetaData");
+
+    if (mesh_property_ip_meta_data.getMeshItemType() !=
+        MeshLib::MeshItemType::IntegrationPoint)
+    {
+        OGS_FATAL("IntegrationPointMetaData array must be field data.");
+    }
+
+    // Find the current integration point data entry and extract the
+    // meta data.
+    auto const ip_meta_data = extractIntegrationPointMetaData(
+        json::parse(mesh_property_ip_meta_data.begin(),
+                    mesh_property_ip_meta_data.end()),
+        name);
+
+    return ip_meta_data;
+}
 }  // namespace ProcessLib

ProcessLib/Output/IntegrationPointWriter.h

@@ -9,6 +9,7 @@
  *
  */
 #include <memory>
+#include <nlohmann/json_fwd.hpp>
 #include <vector>
 
 #pragma once
@@ -51,4 +52,7 @@ struct IntegrationPointMetaData
     int const integration_order;
 };
 
+// Get integration point meta data.
+IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
+                                                     std::string const& name);
 }  // namespace ProcessLib

ProcessLib/SmallDeformation/LocalAssemblerInterface.h

@@ -26,6 +26,10 @@ struct SmallDeformationLocalAssemblerInterface
       public MaterialForcesInterface,
       public NumLib::ExtrapolatableElement
 {
+    virtual std::size_t setIPDataInitialConditions(
+        std::string const& name, double const* values,
+        int const integration_order) = 0;
+
     virtual std::vector<double> const& getIntPtFreeEnergyDensity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -151,6 +151,29 @@ class SmallDeformationLocalAssembler
         }
     }
 
+    /// Returns number of read integration points.
+    std::size_t setIPDataInitialConditions(std::string const& name,
+                                           double const* values,
+                                           int const integration_order) override
+    {
+        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")
+        {
+            return setSigma(values);
+        }
+
+        return 0;
+    }
+
     void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
                   std::vector<double>& /*local_M_data*/,
                   std::vector<double>& /*local_K_data*/,
@@ -314,6 +337,29 @@ class SmallDeformationLocalAssembler
         return cache;
     }
 
+    std::size_t setSigma(double const* values)
+    {
+        auto const kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        auto const n_integration_points = _ip_data.size();
+
+        std::vector<double> ip_sigma_values;
+        auto sigma_values =
+            Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
+                                     Eigen::ColMajor> const>(
+                values, kelvin_vector_size, n_integration_points);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            _ip_data[ip].sigma =
+                MathLib::KelvinVector::symmetricTensorToKelvinVector(
+                    sigma_values.col(ip));
+        }
+
+        return n_integration_points;
+    }
+
     // TODO (naumov) This method is same as getIntPtSigma but for arguments and
     // the ordering of the cache_mat.
     // There should be only one.

ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include <cassert>
+#include <nlohmann/json.hpp>
 
 #include "BaseLib/Functional.h"
 #include "ProcessLib/Process.h"
@@ -72,6 +73,8 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
+    using nlohmann::json;
+
     ProcessLib::SmallDeformation::createLocalAssemblers<
         DisplacementDim, SmallDeformationLocalAssembler>(
         mesh.getElements(), dof_table, _local_assemblers,
@@ -157,6 +160,56 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
             makeExtrapolator(num_components, getExtrapolator(),
                              _local_assemblers, std::move(getIntPtValues)));
     }
+
+    // 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& mesh_property =
+            *mesh.getProperties().template getPropertyVector<double>(name);
+
+        // The mesh property must be defined on integration points.
+        if (mesh_property.getMeshItemType() !=
+            MeshLib::MeshItemType::IntegrationPoint)
+        {
+            continue;
+        }
+
+        auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
+
+        // Check the number of components.
+        if (ip_meta_data.n_components != mesh_property.getNumberOfComponents())
+        {
+            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.c_str(),
+                mesh_property.getNumberOfComponents());
+        }
+
+        // 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, &mesh_property[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;
+        }
+    }
 }
 
 template <int DisplacementDim>