[PL/SD] Use media replacing solid_density

The evaluation of solid density is moved to the
updateConstitutiveRelations.

Documentation/ProjectFile/prj/processes/process/SMALL_DEFORMATION/t_solid_density.md

-\copydoc ProcessLib::SmallDeformation::SmallDeformationProcessData::solid_density
+No longer used as of version 6.5.0. Replaced with media properties for the solid
+phase. See corresponding [merge-request 4739](https://gitlab.opengeosys.org/ogs/ogs/-/merge_requests/4739)
+for details and a conversion script for the project files.

ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp

@@ -12,6 +12,7 @@
 
 #include <cassert>
 
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
@@ -23,6 +24,16 @@ namespace ProcessLib
 {
 namespace SmallDeformation
 {
+void checkMPLProperties(
+    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
+{
+    for (auto const& m : media)
+    {
+        checkRequiredProperties(m.second->phase("Solid"),
+                                {{MaterialPropertyLib::density}});
+    }
+}
+
 template <int DisplacementDim>
 std::unique_ptr<Process> createSmallDeformationProcess(
     std::string const& name,
@@ -73,12 +84,13 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, local_coordinate_system, config);
 
-    // Solid density
-    auto const& solid_density = ParameterLib::findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION__solid_density}
-        "solid_density", parameters, 1, &mesh);
-    DBUG("Use '{:s}' as solid density parameter.", solid_density.name);
+    //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__solid_density}
+    if (config.getConfigParameterOptional<std::string>("solid_density"))
+    {
+        OGS_FATAL(
+            "The <solid_density> tag has been removed. Use <media> definitions "
+            "to specify solid's density.");
+    }
 
     // Specific body force
     Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
@@ -99,6 +111,14 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
 
+    auto media_map =
+        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+    DBUG(
+        "Check the media properties of SmallDeformation process "
+        "...");
+    checkMPLProperties(media);
+    DBUG("Media properties verified.");
+
     // Reference temperature
     auto const reference_temperature = ParameterLib::findOptionalTagParameter<
         double>(
@@ -119,9 +139,12 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         &mesh);
 
     SmallDeformationProcessData<DisplacementDim> process_data{
-        materialIDs(mesh),   std::move(solid_constitutive_relations),
-        initial_stress,      solid_density,
-        specific_body_force, reference_temperature};
+        materialIDs(mesh),
+        std::move(media_map),
+        std::move(solid_constitutive_relations),
+        initial_stress,
+        specific_body_force,
+        reference_temperature};
 
     SecondaryVariableCollection secondary_variables;
 

ProcessLib/SmallDeformation/SmallDeformationFEM.h

@@ -262,9 +262,9 @@ public:
         }
     }
 
-    // Updates sigma, eps, and state through passed integration point data
-    // Returns tangent stiffness.
-    MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
+    /// Updates sigma, eps, and state through passed integration point data.
+    /// \returns tangent stiffness and density.
+    std::tuple<MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>, double>
     updateConstitutiveRelations(
         Eigen::Ref<Eigen::VectorXd const> const& u,
         Eigen::Ref<Eigen::VectorXd const> const& u_prev,
@@ -273,6 +273,10 @@ public:
         IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>&
             ip_data) const
     {
+        auto const& solid_phase =
+            this->_process_data.media_map->getMedium(this->_element.getID())
+                ->phase("Solid");
+
         MPL::VariableArray variables_prev;
         MPL::VariableArray variables;
 
@@ -313,6 +317,10 @@ public:
                 eps);
         variables.temperature = T_ref;
 
+        auto const rho =
+            solid_phase[MPL::PropertyType::density].template value<double>(
+                variables, x_position, t, dt);
+
         auto&& solution = ip_data.solid_material.integrateStress(
             variables_prev, variables, t, x_position, dt, *state);
 
@@ -324,7 +332,7 @@ public:
         MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
         std::tie(sigma, state, C) = std::move(*solution);
 
-        return C;
+        return {C, rho};
     }
 
     void assemble(double const /*t*/, double const /*dt*/,
@@ -383,10 +391,9 @@ public:
 
             auto const& sigma = _ip_data[ip].sigma;
 
-            auto const C = updateConstitutiveRelations(u, u_prev, x_position, t,
-                                                       dt, _ip_data[ip]);
+            auto const [C, rho] = updateConstitutiveRelations(
+                u, u_prev, x_position, t, dt, _ip_data[ip]);
 
-            auto const rho = _process_data.solid_density(t, x_position)[0];
             local_b.noalias() -=
                 (B.transpose() * sigma - N_u_op(N).transpose() * rho * b) * w;
             local_Jac.noalias() += B.transpose() * C * B * w;

ProcessLib/SmallDeformation/SmallDeformationProcessData.h

@@ -14,6 +14,7 @@
 #include <memory>
 #include <utility>
 
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "ParameterLib/Parameter.h"
 
 namespace MaterialLib
@@ -33,6 +34,9 @@ struct SmallDeformationProcessData
 {
     MeshLib::PropertyVector<int> const* const material_ids = nullptr;
 
+    std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+        media_map = nullptr;
+
     std::map<
         int,
         std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
@@ -42,9 +46,6 @@ struct SmallDeformationProcessData
     /// representation of length 4 or 6, ParameterLib::Parameter<double>.
     ParameterLib::Parameter<double> const* const initial_stress;
 
-    /// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
-    ParameterLib::Parameter<double> const& solid_density;
-
     /// Specific body forces applied to the solid.
     /// It is usually used to apply gravitational forces.
     /// A vector of displacement dimension's length.