Merge branch 'u0_function_arg' into 'master'

Pass the solution vectors of all coupled processes to setInitialConditions

See merge request ogs/ogs!4950

NumLib/DOF/DOFTableUtil.cpp

@@ -10,8 +10,11 @@
 
 #include "DOFTableUtil.h"
 
+#include <algorithm>
 #include <cassert>
+#include <functional>
 
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 namespace NumLib
 {
 namespace
@@ -177,4 +180,39 @@ double norm(GlobalVector const& x, unsigned const global_component,
     }
 }
 
+std::vector<NumLib::LocalToGlobalIndexMap const*> getDOFTables(
+    int const number_of_processes,
+    std::function<NumLib::LocalToGlobalIndexMap const&(const int)>
+        get_single_dof_table)
+{
+    std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+    dof_tables.reserve(number_of_processes);
+    std::generate_n(std::back_inserter(dof_tables), number_of_processes,
+                    [&]() { return &get_single_dof_table(dof_tables.size()); });
+    return dof_tables;
+}
+
+Eigen::VectorXd getLocalX(
+    std::size_t const mesh_item_id,
+    std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
+    std::vector<GlobalVector*> const& x)
+{
+    Eigen::VectorXd local_x_vec;
+
+    auto const n_processes = x.size();
+    for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
+    {
+        auto const indices =
+            NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
+        assert(!indices.empty());
+        auto const last = local_x_vec.size();
+        local_x_vec.conservativeResize(last + indices.size());
+        auto const local_solution = x[process_id]->get(indices);
+        assert(indices.size() == local_solution.size());
+        local_x_vec.tail(local_solution.size()).noalias() =
+            MathLib::toVector(local_solution);
+    }
+    return local_x_vec;
+}
+
 }  // namespace NumLib

NumLib/DOF/DOFTableUtil.h

@@ -10,6 +10,8 @@
 
 #pragma once
 
+#include <Eigen/Core>
+
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/LinAlgEnums.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
@@ -143,4 +145,14 @@ void transformVariableFromGlobalVector(
         }
     }
 }
+
+std::vector<NumLib::LocalToGlobalIndexMap const*> getDOFTables(
+    int const number_of_processes,
+    std::function<NumLib::LocalToGlobalIndexMap const&(const int)>
+        get_single_dof_table);
+
+Eigen::VectorXd getLocalX(
+    std::size_t const mesh_item_id,
+    std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
+    std::vector<GlobalVector*> const& x);
 }  // namespace NumLib

ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h

@@ -815,78 +815,57 @@ template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
           int DisplacementDim>
 void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                   ShapeFunctionPressure, DisplacementDim>::
-    setInitialConditionsConcrete(std::vector<double> const& local_x,
+    setInitialConditionsConcrete(Eigen::VectorXd const local_x,
                                  double const t,
-                                 bool const use_monolithic_scheme,
-                                 int const process_id)
+                                 int const /*process_id*/)
 {
-    if (!use_monolithic_scheme &&
-        process_id == _process_data.hydraulic_process_id)
-    {
-        return;
-    }
+    ParameterLib::SpatialPosition x_position;
+    x_position.setElementID(_element.getID());
 
-    if (use_monolithic_scheme ||
-        process_id == _process_data.mechanics_related_process_id)
-    {
-        ParameterLib::SpatialPosition x_position;
-        x_position.setElementID(_element.getID());
+    auto const& medium = _process_data.media_map.getMedium(_element.getID());
 
-        auto const& medium =
-            _process_data.media_map.getMedium(_element.getID());
+    auto const p = local_x.template segment<pressure_size>(pressure_index);
+    auto const u =
+        local_x.template segment<displacement_size>(displacement_index);
 
-        int const displacement_offset =
-            use_monolithic_scheme ? displacement_index : 0;
+    auto const& identity2 = Invariants::identity2;
+    const double dt = 0.0;
 
-        auto const u =
-            Eigen::Map<typename ShapeMatricesTypeDisplacement::
-                           template VectorType<displacement_size> const>(
-                local_x.data() + displacement_offset, displacement_size);
+    MPL::VariableArray vars;
 
-        int const pressure_offset = use_monolithic_scheme ? pressure_index : 0;
+    int const n_integration_points = _integration_method.getNumberOfPoints();
+    for (int ip = 0; ip < n_integration_points; ip++)
+    {
+        x_position.setIntegrationPoint(ip);
+        auto const& N_u = _ip_data[ip].N_u;
+        auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
-        auto const p =
-            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-                pressure_size> const>(local_x.data() + pressure_offset,
-                                      pressure_size);
-        auto const& identity2 = Invariants::identity2;
-        const double dt = 0.0;
+        auto const x_coord =
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
+        auto const B =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                dNdx_u, N_u, x_coord, _is_axially_symmetric);
 
-        MPL::VariableArray vars;
+        auto& eps = _ip_data[ip].eps;
+        eps.noalias() = B * u;
+        vars.mechanical_strain
+            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+                eps);
 
-        int const n_integration_points =
-            _integration_method.getNumberOfPoints();
-        for (int ip = 0; ip < n_integration_points; ip++)
+        if (_process_data.initial_stress.isTotalStress())
         {
-            x_position.setIntegrationPoint(ip);
-            auto const& N_u = _ip_data[ip].N_u;
-            auto const& dNdx_u = _ip_data[ip].dNdx_u;
-
-            auto const x_coord =
-                NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
-                                               ShapeMatricesTypeDisplacement>(
-                    _element, N_u);
-            auto const B = LinearBMatrix::computeBMatrix<
-                DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
-                typename BMatricesType::BMatrixType>(dNdx_u, N_u, x_coord,
-                                                     _is_axially_symmetric);
-
-            auto& eps = _ip_data[ip].eps;
-            eps.noalias() = B * u;
-            vars.mechanical_strain.emplace<
-                MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(eps);
-
-            if (_process_data.initial_stress.isTotalStress())
-            {
-                auto const& N_p = _ip_data[ip].N_p;
-                auto const alpha_b =
-                    medium->property(MPL::PropertyType::biot_coefficient)
-                        .template value<double>(vars, x_position, t, dt);
+            auto const& N_p = _ip_data[ip].N_p;
+            auto const alpha_b =
+                medium->property(MPL::PropertyType::biot_coefficient)
+                    .template value<double>(vars, x_position, t, dt);
 
-                auto& sigma_eff = _ip_data[ip].sigma_eff;
-                sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
-                _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
-            }
+            auto& sigma_eff = _ip_data[ip].sigma_eff;
+            sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
+            _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
         }
     }
 }

ProcessLib/HydroMechanics/HydroMechanicsFEM.h

@@ -276,9 +276,8 @@ class HydroMechanicsLocalAssembler
                                      double const t, double const dt,
                                      int const process_id) override;
 
-    void setInitialConditionsConcrete(std::vector<double> const& local_x,
+    void setInitialConditionsConcrete(Eigen::VectorXd const local_x,
                                       double const t,
-                                      bool const use_monolithic_scheme,
                                       int const process_id) override;
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(

ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp

@@ -17,6 +17,7 @@
 #include "MeshLib/Elements/Utils.h"
 #include "MeshLib/Utils/getOrCreateMeshProperty.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
@@ -398,17 +399,16 @@ void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
     }
 
     DBUG("PostTimestep HydroMechanicsProcess.");
-    std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
-    auto const n_processes = x.size();
-    dof_tables.reserve(n_processes);
-    for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
-    {
-        dof_tables.push_back(&getDOFTable(process_id));
-    }
 
+    auto get_a_dof_table_func = [this](const int process_id) -> auto&
+    {
+        return getDOFTable(process_id);
+    };
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerIF::postTimestep, _local_assemblers,
-        getActiveElementIDs(), dof_tables, x, x_prev, t, dt, process_id);
+        getActiveElementIDs(),
+        NumLib::getDOFTables(x.size(), get_a_dof_table_func), x, x_prev, t, dt,
+        process_id);
 }
 
 template <int DisplacementDim>
@@ -432,17 +432,22 @@ void HydroMechanicsProcess<DisplacementDim>::
                                         double const t,
                                         int const process_id)
 {
-    if (process_id != _process_data.hydraulic_process_id)
+    // So far, this function only sets the initial stress using the input data.
+    if (process_id != _process_data.mechanics_related_process_id)
     {
         return;
     }
 
     DBUG("Set initial conditions of HydroMechanicsProcess.");
 
+    auto get_a_dof_table_func = [this](const int process_id) -> auto&
+    {
+        return getDOFTable(process_id);
+    };
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerIF::setInitialConditions, _local_assemblers,
-        getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
-        _process_data.isMonolithicSchemeUsed(), process_id);
+        getActiveElementIDs(),
+        NumLib::getDOFTables(x.size(), get_a_dof_table_func), x, t, process_id);
 }
 
 template <int DisplacementDim>
@@ -456,17 +461,16 @@ void HydroMechanicsProcess<DisplacementDim>::computeSecondaryVariableConcrete(
     }
 
     DBUG("Compute the secondary variables for HydroMechanicsProcess.");
-    std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
-    auto const n_processes = x.size();
-    dof_tables.reserve(n_processes);
-    for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
-    {
-        dof_tables.push_back(&getDOFTable(process_id));
-    }
 
+    auto get_a_dof_table_func = [this](const int process_id) -> auto&
+    {
+        return getDOFTable(process_id);
+    };
     GlobalExecutor::executeSelectedMemberOnDereferenced(
         &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
-        getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id);
+        getActiveElementIDs(),
+        NumLib::getDOFTables(x.size(), get_a_dof_table_func), t, dt, x, x_prev,
+        process_id);
 }
 
 template <int DisplacementDim>

ProcessLib/HydroMechanics/Tests.cmake

@@ -770,6 +770,22 @@ AddTest(
     total_initial_stress_HM_ts_1_t_1000000.000000.vtu total_initial_stress_HM_ts_1_t_1000000.000000.vtu sigma sigma 1.0e-8 1.e-10
 )
 
+AddTest(
+    NAME HydroMechanicsTotalInitialStress_Staggered
+    PATH HydroMechanics/TotalInitialStress
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS total_initial_stress_HM_staggered.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    RUNTIME 1
+    DIFF_DATA
+    total_initial_stress_HM_ts_0_t_0.000000.vtu total_initial_stress_HM_staggered_ts_0_t_0.000000.vtu pressure pressure 1.0e-10 1.e-10
+    total_initial_stress_HM_ts_0_t_0.000000.vtu total_initial_stress_HM_staggered_ts_0_t_0.000000.vtu sigma sigma 1.0e-8 1.e-10
+    total_initial_stress_HM_ts_1_t_1000000.000000.vtu total_initial_stress_HM_staggered_ts_1_t_1000000.000000.vtu pressure pressure 1.0e-10 1.e-10
+    total_initial_stress_HM_ts_1_t_1000000.000000.vtu total_initial_stress_HM_staggered_ts_1_t_1000000.000000.vtu sigma sigma 1.0e-8 1.e-10
+)
+
 #Parallel computing with PETSc enabled.
 AddTest(
     NAME ParallelFEM_SimpleHM_SingleMesh

ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp

@@ -378,20 +378,18 @@ void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
     if (process_id == 0)
     {
         DBUG("PostTimestep HydroMechanicsProcess.");
-        std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
-        auto const n_processes = x.size();
-        dof_tables.reserve(n_processes);
-        for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
-        {
-            dof_tables.push_back(&getDOFTable(process_id));
-        }
 
+        auto get_a_dof_table_func = [this](const int processe_id) -> auto&
+        {
+            return getDOFTable(processe_id);
+        };
         ProcessLib::ProcessVariable const& pv =
             getProcessVariables(process_id)[0];
         GlobalExecutor::executeSelectedMemberOnDereferenced(
             &HydroMechanicsLocalAssemblerInterface::postTimestep,
-            _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, x_prev,
-            t, dt, process_id);
+            _local_assemblers, pv.getActiveElementIDs(),
+            NumLib::getDOFTables(x.size(), get_a_dof_table_func), x, x_prev, t,
+            dt, process_id);
     }
 
     DBUG("Compute the secondary variables for HydroMechanicsProcess.");