[PL] SD: Set integration point initial conditions.

ProcessLib/SmallDeformation/LocalAssemblerInterface.h

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

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -151,6 +151,15 @@ class SmallDeformationLocalAssembler
         }
     }
 
+    void setIPDataInitialConditions(std::string const& name,
+                                    double const* values) override
+    {
+        if (name == "sigma_ip")
+        {
+            setSigma(values);
+        }
+    }
+
     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 +323,27 @@ class SmallDeformationLocalAssembler
         return cache;
     }
 
+    void 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));
+        }
+    }
+
     std::vector<double> getSigma() const override
     {
         using KelvinVectorType = typename BMatricesType::KelvinVectorType;

ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h

@@ -157,6 +157,60 @@ 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)
+    {
+        auto const& name = ip_writer->name();
+        if (!mesh.getProperties().existsPropertyVector<double>(name))
+        {
+            continue;
+        }
+
+        auto const& mesh_property =
+            *mesh.getProperties().template getPropertyVector<double>(name);
+
+        if (mesh_property.getMeshItemType() !=
+            MeshLib::MeshItemType::IntegrationPoint)
+        {
+            continue;
+        }
+
+        auto const offsets_array_name = name + "_offsets";
+        if (!mesh.getProperties().existsPropertyVector<std::size_t>(
+                offsets_array_name))
+        {
+            OGS_FATAL(
+                "Integration point data '%s' is present in the vtk field data "
+                "but the corresponding '%s' array is not available.",
+                name.c_str(), offsets_array_name.c_str());
+        }
+
+        auto const& mesh_property_offsets =
+            *mesh.getProperties().template getPropertyVector<std::size_t>(
+                offsets_array_name);
+
+        if (mesh_property_offsets.getMeshItemType() !=
+            MeshLib::MeshItemType::Cell)
+        {
+            OGS_FATAL(
+                "Integration point data '%s' is present in the vtk field data "
+                "but the corresponding '%s' array is not defined on cells.",
+                name.c_str(), offsets_array_name.c_str());
+        }
+
+        // Now we have a properly named vtk's field data array and the
+        // corresponding offsets array.
+        for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+        {
+            auto& local_asm = _local_assemblers[i];
+            auto const offset = mesh_property_offsets[i];
+
+            // TODO (naumov) Check sizes / read size / etc.
+            // OR reconstruct dimensions from size / component = ip_points
+            local_asm->setIPDataInitialConditions(name, &mesh_property[offset]);
+        }
+    }
 }
 
 template <int DisplacementDim>