[PL] Implementation of RichardsMechanics process.

Applications/ApplicationsLib/ProjectData.cpp

@@ -72,6 +72,9 @@
 #ifdef OGS_BUILD_PROCESS_RICHARDSFLOW
 #include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
 #endif
+#ifdef OGS_BUILD_PROCESS_RICHARDSMECHANICS
+#include "ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h"
+#endif
 #ifdef OGS_BUILD_PROCESS_SMALLDEFORMATION
 #include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
 #endif
@@ -554,6 +557,30 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
         }
         else
 #endif
+#ifdef OGS_BUILD_PROCESS_RICHARDSMECHANICS
+            if (type == "RICHARDS_MECHANICS")
+        {
+            //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__dimension}
+            switch (process_config.getConfigParameter<int>("dimension"))
+            {
+                case 2:
+                    process = ProcessLib::RichardsMechanics::
+                        createRichardsMechanicsProcess<2>(
+                            *_mesh_vec[0], std::move(jacobian_assembler),
+                            _process_variables, _parameters, integration_order,
+                            process_config);
+                    break;
+                case 3:
+                    process = ProcessLib::RichardsMechanics::
+                        createRichardsMechanicsProcess<3>(
+                            *_mesh_vec[0], std::move(jacobian_assembler),
+                            _process_variables, _parameters, integration_order,
+                            process_config);
+                    break;
+            }
+        }
+        else
+#endif
 #ifdef OGS_BUILD_PROCESS_TWOPHASEFLOWWITHPP
             if (type == "TWOPHASE_FLOW_PP")
         {

CMakeLists.txt

@@ -99,6 +99,7 @@ set(ProcessesList
     PhaseField
     RichardsComponentTransport
     RichardsFlow
+    RichardsMechanics
     SmallDeformation
     TES
     ThermalTwoPhaseFlowWithPP

ProcessLib/RichardsMechanics/CMakeLists.txt

+APPEND_SOURCE_FILES(SOURCES)
+
+add_library(RichardsMechanics ${SOURCES})
+target_link_libraries(RichardsMechanics PUBLIC ProcessLib RichardsFlow)
+
+include(Tests.cmake)

ProcessLib/RichardsMechanics/CreateLocalAssemblers.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <vector>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "LocalDataInitializer.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+namespace detail
+{
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    // Shape matrices initializer
+    using LocalDataInitializer =
+        LocalDataInitializer<LocalAssemblerInterface,
+                             LocalAssemblerImplementation, GlobalDim,
+                             ExtraCtorArgs...>;
+
+    DBUG("Create local assemblers.");
+    // Populate the vector of local assemblers.
+    local_assemblers.resize(mesh_elements.size());
+
+    LocalDataInitializer initializer(dof_table, shapefunction_order);
+
+    DBUG("Calling local assembler builder for all mesh elements.");
+    GlobalExecutor::transformDereferenced(
+        initializer, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace detail
+
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ *         Those arguments will be passed to the constructor of
+ *         \c LocalAssemblerImplementation.
+ *
+ * The first two template parameters cannot be deduced from the arguments.
+ * Therefore they always have to be provided manually.
+ */
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    const unsigned /*dimension*/,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    DBUG("Create local assemblers.");
+
+    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
+        dof_table, shapefunction_order, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+}  // namespace RichardsMechanics
+
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateRichardsMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "RichardsMechanicsProcess.h"
+#include "RichardsMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createRichardsMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{prj__processes__process__type}
+    config.checkConfigParameter("type", "RICHARDS_MECHANICS");
+    DBUG("Create RichardsMechanicsProcess.");
+
+    auto const staggered_scheme =
+        //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__coupling_scheme}
+        config.getConfigParameterOptional<std::string>("coupling_scheme");
+    const bool use_monolithic_scheme =
+        !(staggered_scheme && (*staggered_scheme == "staggered"));
+
+    // Process variable.
+
+    //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__process_variables}
+    auto const pv_config = config.getConfigSubtree("process_variables");
+
+    ProcessVariable* variable_p;
+    ProcessVariable* variable_u;
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+        process_variables;
+    if (use_monolithic_scheme)  // monolithic scheme.
+    {
+        auto per_process_variables = findProcessVariables(
+            variables, pv_config,
+            {//! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__process_variables__pressure}
+             "pressure",
+             //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__process_variables__displacement}
+             "displacement"});
+        variable_p = &per_process_variables[0].get();
+        variable_u = &per_process_variables[1].get();
+        process_variables.push_back(std::move(per_process_variables));
+    }
+    else  // staggered scheme.
+    {
+        using namespace std::string_literals;
+        for (auto const& variable_name : {"pressure"s, "displacement"s})
+        {
+            auto per_process_variables =
+                findProcessVariables(variables, pv_config, {variable_name});
+            process_variables.push_back(std::move(per_process_variables));
+        }
+        variable_p = &process_variables[0][0].get();
+        variable_u = &process_variables[1][0].get();
+    }
+
+    DBUG("Associate displacement with process variable \'%s\'.",
+         variable_u->getName().c_str());
+
+    if (variable_u->getNumberOfComponents() != DisplacementDim)
+    {
+        OGS_FATAL(
+            "Number of components of the process variable '%s' is different "
+            "from the displacement dimension: got %d, expected %d",
+            variable_u->getName().c_str(),
+            variable_u->getNumberOfComponents(),
+            DisplacementDim);
+    }
+
+    DBUG("Associate pressure with process variable \'%s\'.",
+         variable_p->getName().c_str());
+    if (variable_p->getNumberOfComponents() != 1)
+    {
+        OGS_FATAL(
+            "Pressure process variable '%s' is not a scalar variable but has "
+            "%d components.",
+            variable_p->getName().c_str(),
+            variable_p->getNumberOfComponents());
+    }
+
+    // Constitutive relation.
+    auto solid_material =
+        MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+            parameters, config);
+
+    // Intrinsic permeability
+    auto& intrinsic_permeability = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__intrinsic_permeability}
+        "intrinsic_permeability", parameters, 1);
+
+    DBUG("Use \'%s\' as intrinsic conductivity parameter.",
+         intrinsic_permeability.name.c_str());
+
+    // Fluid bulk modulus
+    auto& fluid_bulk_modulus = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__fluid_bulk_modulus}
+        "fluid_bulk_modulus", parameters, 1);
+    DBUG("Use \'%s\' as fluid bulk modulus parameter.",
+         fluid_bulk_modulus.name.c_str());
+
+    // Biot coefficient
+    auto& biot_coefficient = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__biot_coefficient}
+        "biot_coefficient", parameters, 1);
+    DBUG("Use \'%s\' as Biot coefficient parameter.",
+         biot_coefficient.name.c_str());
+
+    // Solid density
+    auto& solid_density = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__solid_density}
+        "solid_density", parameters, 1);
+    DBUG("Use \'%s\' as solid density parameter.", solid_density.name.c_str());
+
+    // Solid bulk modulus
+    auto& solid_bulk_modulus = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__solid_bulk_modulus}
+        "solid_bulk_modulus", parameters, 1);
+    DBUG("Use \'%s\' as solid bulk modulus parameter.",
+         solid_bulk_modulus.name.c_str());
+
+    // Specific body force
+    Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+    {
+        std::vector<double> const b =
+            //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__specific_body_force}
+            config.getConfigParameter<std::vector<double>>(
+                "specific_body_force");
+        if (b.size() != DisplacementDim)
+        {
+            OGS_FATAL(
+                "The size of the specific body force vector does not match the "
+                "displacement dimension. Vector size is %d, displacement "
+                "dimension is %d",
+                b.size(), DisplacementDim);
+        }
+
+        std::copy_n(b.data(), b.size(), specific_body_force.data());
+    }
+
+    auto& temperature = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__temperature}
+        "temperature", parameters, 1);
+
+    auto const& flow_material_config =
+        //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__material_property}
+        config.getConfigSubtree("material_property");
+
+    boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
+    if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+    {
+        INFO(
+            "MaterialIDs vector found; the Richards flow is in heterogeneous "
+            "porous media.");
+        material_ids =
+            *mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+    }
+    else
+    {
+        INFO(
+            "MaterialIDs vector not found; the Richards flow is in homogeneous "
+            "porous media.");
+    }
+    auto flow_material =
+        ProcessLib::RichardsFlow::createRichardsFlowMaterialProperties(
+            flow_material_config, material_ids, parameters);
+
+    RichardsMechanicsProcessData<DisplacementDim> process_data{
+        std::move(flow_material),
+        std::move(solid_material),
+        intrinsic_permeability,
+        fluid_bulk_modulus,
+        biot_coefficient,
+        solid_density,
+        solid_bulk_modulus,
+        temperature,
+        specific_body_force};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"RichardsMechanics_displacement"});
+
+    ProcessLib::createSecondaryVariables(config, secondary_variables,
+                                         named_function_caller);
+
+    return std::make_unique<RichardsMechanicsProcess<DisplacementDim>>(
+        mesh, std::move(jacobian_assembler), parameters, integration_order,
+        std::move(process_variables), std::move(process_data),
+        std::move(secondary_variables), std::move(named_function_caller),
+        use_monolithic_scheme);
+}
+
+template std::unique_ptr<Process> createRichardsMechanicsProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<Process> createRichardsMechanicsProcess<3>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <map>
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace MathLib
+{
+class PiecewiseLinearInterpolation;
+}
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+struct ParameterBase;
+class Process;
+class ProcessVariable;
+}  // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createRichardsMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<Process> createRichardsMechanicsProcess<2>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<Process> createRichardsMechanicsProcess<3>(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<ProcessVariable> const& variables,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    BaseLib::ConfigTree const& config);
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/IntegrationPointData.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <typename BMatricesType, typename ShapeMatrixTypeDisplacement,
+          typename ShapeMatricesTypePressure, int DisplacementDim, int NPoints>
+struct IntegrationPointData final
+{
+    explicit IntegrationPointData(
+        MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material)
+        : solid_material(solid_material),
+          material_state_variables(
+              solid_material.createMaterialStateVariables())
+    {
+        // Initialize current time step values
+        static const int kelvin_vector_size =
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value;
+        sigma_eff.setZero(kelvin_vector_size);
+        eps.setZero(kelvin_vector_size);
+
+        // Previous time step values are not initialized and are set later.
+        eps_prev.resize(kelvin_vector_size);
+        sigma_eff_prev.resize(kelvin_vector_size);
+    }
+
+    typename ShapeMatrixTypeDisplacement::template MatrixType<
+        DisplacementDim, NPoints * DisplacementDim>
+        N_u_op;
+    typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
+    typename BMatricesType::KelvinVectorType eps, eps_prev;
+
+    typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+    typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
+
+    typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
+    typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
+
+    double saturation;
+
+    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables>
+        material_state_variables;
+    double integration_weight;
+
+    void pushBackState()
+    {
+        eps_prev = eps;
+        sigma_eff_prev = sigma_eff;
+        material_state_variables->pushBackState();
+    }
+
+    template <typename DisplacementVectorType>
+    typename BMatricesType::KelvinMatrixType updateConstitutiveRelation(
+        double const t,
+        SpatialPosition const& x_position,
+        double const dt,
+        DisplacementVectorType const& /*u*/,
+        double const temperature)
+    {
+        auto&& solution = solid_material.integrateStress(
+            t, x_position, dt, eps_prev, eps, sigma_eff_prev,
+            *material_state_variables, temperature);
+
+        if (!solution)
+            OGS_FATAL("Computation of local constitutive relation failed.");
+
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
+        std::tie(sigma_eff, material_state_variables, C) = std::move(*solution);
+
+        return C;
+    }
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/LocalAssemblerInterface.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
+                                 public NumLib::ExtrapolatableElement
+{
+    virtual std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+
+    virtual std::vector<double> const& getIntPtSaturation(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const = 0;
+};
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/LocalDataInitializer.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/FiniteElement/LowerDimShapeTable.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+
+#ifndef OGS_MAX_ELEMENT_DIM
+static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
+#endif
+
+#ifndef OGS_MAX_ELEMENT_ORDER
+static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
+#endif
+
+// The following macros decide which element types will be compiled, i.e.
+// which element types will be available for use in simulations.
+
+#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
+#else
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_CUBOID
+#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
+#else
+#define ENABLED_ELEMENT_TYPE_CUBOID 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PRISM
+#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
+#else
+#define ENABLED_ELEMENT_TYPE_PRISM 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PYRAMID
+#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
+#else
+#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
+#endif
+
+// Dependent element types.
+// Faces of tets, pyramids and prisms are triangles
+#define ENABLED_ELEMENT_TYPE_TRI                                       \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of hexes, pyramids and prisms are quads
+#define ENABLED_ELEMENT_TYPE_QUAD                                     \
+    ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+
+// All enabled element types
+#define OGS_ENABLED_ELEMENTS                                          \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
+     (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
+
+// Include only what is needed (Well, the conditions are not sharp).
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
+#endif
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+/// The LocalDataInitializer is a functor creating a local assembler data with
+/// corresponding to the mesh element type shape functions and calling
+/// initialization of the new local assembler data.
+/// For example for MeshLib::Quad a local assembler data with template argument
+/// NumLib::ShapeQuad4 is created.
+///
+/// \attention This is modified version of the ProcessLib::LocalDataInitializer
+/// class which does not include line elements, allows only shapefunction of
+/// order 2.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, int>
+          class RichardsMechanicsLocalAssembler,
+          int GlobalDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+    LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
+                         const unsigned shapefunction_order)
+        : _dof_table(dof_table)
+    {
+        if (shapefunction_order != 2)
+            OGS_FATAL(
+                "The given shape function order %d is not supported.\nOnly "
+                "shape functions of order 2 are supported.",
+                shapefunction_order);
+            // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Quad8))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+        _builder[std::type_index(typeid(MeshLib::Quad9))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Hex20))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+        // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Tri6))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Tet10))] =
+            makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+        // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Prism15))] =
+            makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+        // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+        _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+            makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
+#endif
+    }
+
+    /// Sets the provided \c data_ptr to the newly created local assembler data.
+    ///
+    /// \attention
+    /// The index \c id is not necessarily the mesh item's id. Especially when
+    /// having multiple meshes it will differ from the latter.
+    void operator()(std::size_t const id,
+                    MeshLib::Element const& mesh_item,
+                    LADataIntfPtr& data_ptr,
+                    ConstructorArgs&&... args) const
+    {
+        auto const type_idx = std::type_index(typeid(mesh_item));
+        auto const it = _builder.find(type_idx);
+
+        if (it != _builder.end())
+        {
+            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
+            data_ptr = it->second(mesh_item, num_local_dof,
+                                  std::forward<ConstructorArgs>(args)...);
+        }
+        else
+        {
+            OGS_FATAL(
+                "You are trying to build a local assembler for an unknown mesh "
+                "element type (%s)."
+                " Maybe you have disabled this mesh element type in your build "
+                "configuration or this process requires higher order elements.",
+                type_idx.name());
+        }
+    }
+
+private:
+    using LADataBuilder =
+        std::function<LADataIntfPtr(MeshLib::Element const& e,
+                                    std::size_t const local_matrix_size,
+                                    ConstructorArgs&&...)>;
+
+    template <typename ShapeFunctionDisplacement>
+    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+        typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
+
+    template <typename ShapeFunctionDisplacement,
+              typename ShapeFunctionPressure>
+    using LAData = RichardsMechanicsLocalAssembler<
+        ShapeFunctionDisplacement, ShapeFunctionPressure,
+        IntegrationMethod<ShapeFunctionDisplacement>, GlobalDim>;
+
+    /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+    /// depending whether the global dimension is less than the shape function's
+    /// dimension or not.
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder()
+    {
+        return makeLocalAssemblerBuilder<ShapeFunctionDisplacement>(
+            static_cast<std::integral_constant<
+                bool, (GlobalDim >= ShapeFunctionDisplacement::DIM)>*>(
+                nullptr));
+    }
+
+    /// Mapping of element types to local assembler constructors.
+    std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+    // local assembler builder implementations.
+private:
+    /// Generates a function that creates a new LocalAssembler of type
+    /// LAData<ShapeFunctionDisplacement>. Only functions with shape function's
+    /// dimension less or equal to the global dimension are instantiated, e.g.
+    /// following combinations of shape functions and global dimensions: (Line2,
+    /// 1),
+    /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
+    {
+        // (Lower order elements = Order(ShapeFunctionDisplacement) - 1).
+        using ShapeFunctionPressure =
+            typename NumLib::LowerDim<ShapeFunctionDisplacement>::type;
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  ConstructorArgs&&... args) {
+            return LADataIntfPtr{
+                new LAData<ShapeFunctionDisplacement, ShapeFunctionPressure>{
+                    e, local_matrix_size,
+                    std::forward<ConstructorArgs>(args)...}};
+        };
+    }
+
+    /// Returns nullptr for shape functions whose dimensions are less than the
+    /// global dimension.
+    template <typename ShapeFunctionDisplacement>
+    static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
+    {
+        return nullptr;
+    }
+};
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib
+
+#undef ENABLED_ELEMENT_TYPE_SIMPLEX
+#undef ENABLED_ELEMENT_TYPE_CUBOID
+#undef ENABLED_ELEMENT_TYPE_PYRAMID
+#undef ENABLED_ELEMENT_TYPE_PRISM
+#undef ENABLED_ELEMENT_TYPE_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS

ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/KelvinVector.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "RichardsMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+/// Used by for extrapolation of the integration point values. It is ordered
+/// (and stored) by integration points.
+template <typename ShapeMatrixType>
+struct SecondaryData
+{
+    std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N_u;
+};
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+class RichardsMechanicsLocalAssembler : public LocalAssemblerInterface
+{
+public:
+    using ShapeMatricesTypeDisplacement =
+        ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+    using ShapeMatricesTypePressure =
+        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
+
+    static int const KelvinVectorSize =
+        MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+    using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+    RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler const&) =
+        delete;
+    RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler&&) = delete;
+
+    RichardsMechanicsLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const /*local_matrix_size*/,
+        bool const is_axially_symmetric,
+        unsigned const integration_order,
+        RichardsMechanicsProcessData<DisplacementDim>& process_data);
+
+    void assemble(double const t, std::vector<double> const& local_x,
+                  std::vector<double>& local_M_data,
+                  std::vector<double>& local_K_data,
+                  std::vector<double>& local_rhs_data) override;
+
+    void assembleWithJacobian(double const t,
+                              std::vector<double> const& local_x,
+                              std::vector<double> const& local_xdot,
+                              const double /*dxdot_dx*/, const double /*dx_dx*/,
+                              std::vector<double>& /*local_M_data*/,
+                              std::vector<double>& /*local_K_data*/,
+                              std::vector<double>& local_rhs_data,
+                              std::vector<double>& local_Jac_data) override;
+
+    void assembleWithJacobianForStaggeredScheme(
+        double const t, std::vector<double> const& local_xdot,
+        const double dxdot_dx, const double dx_dx,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        LocalCoupledSolutions const& local_coupled_solutions) override;
+
+    void preTimestepConcrete(std::vector<double> const& /*local_x*/,
+                             double const /*t*/,
+                             double const /*delta_t*/) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data[ip].pushBackState();
+        }
+    }
+
+    void computeSecondaryVariableConcrete(
+        double const t, std::vector<double> const& local_x) override;
+
+    void postNonLinearSolverConcrete(std::vector<double> const& local_x,
+                                     double const t,
+                                     bool const use_monolithic_scheme) override;
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N_u = _secondary_data.N_u[integration_point];
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
+    }
+
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override;
+
+    std::vector<double> const& getIntPtSaturation(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& /*cache*/) const override;
+
+    std::vector<double> const& getIntPtSigma(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override;
+
+    std::vector<double> const& getIntPtEpsilon(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& cache) const override;
+
+private:
+    /**
+     * Assemble local matrices and vectors arise from the linearized discretized
+     * weak form of the residual of the momentum balance equation,
+     *      \f[
+     *            \nabla (\sigma - \alpha_b p \mathrm{I}) = f
+     *      \f]
+     * where \f$ \sigma\f$ is the effective stress tensor, \f$p\f$ is the pore
+     * pressure, \f$\alpha_b\f$ is the Biot constant, \f$\mathrm{I}\f$ is the
+     * identity tensor, and \f$f\f$ is the body force.
+     *
+     * @param t               Time
+     * @param local_xdot      Nodal values of \f$\dot{x}\f$ of an element.
+     * @param dxdot_dx        Value of \f$\dot{x} \cdot dx\f$.
+     * @param dx_dx           Value of \f$ x \cdot dx\f$.
+     * @param local_M_data    Mass matrix of an element, which takes the form of
+     *                        \f$ \inta N^T N\mathrm{d}\Omega\f$. Not used.
+     * @param local_K_data    Laplacian matrix of an element, which takes the
+     *         form of \f$ \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\f$.
+     *                        Not used.
+     * @param local_b_data    Right hand side vector of an element.
+     * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
+     *                        method.
+     * @param local_coupled_solutions Nodal values of solutions of the coupled
+     *                                processes of an element.
+     */
+    void assembleWithJacobianForDeformationEquations(
+        double const t, std::vector<double> const& local_xdot,
+        const double dxdot_dx, const double dx_dx,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        LocalCoupledSolutions const& local_coupled_solutions);
+
+    /**
+     * Assemble local matrices and vectors arise from the linearized discretized
+     * weak form of the residual of the mass balance equation of single phase
+     * flow,
+     *      \f[
+     *          \alpha \cdot{p} - \nabla (K (\nabla p + \rho g \nabla z) +
+     *          \alpha_b \nabla \cdot \dot{u}  = Q
+     *      \f]
+     * where \f$ alpha\f$ is a coefficient may depend on storage or the fluid
+     * density change, \f$ \rho\f$ is the fluid density, \f$\g\f$ is the
+     * gravitational acceleration, \f$z\f$ is the vertical coordinate, \f$u\f$
+     * is the displacement, and \f$Q\f$ is the source/sink term.
+     *
+     * @param t               Time
+     * @param local_xdot      Nodal values of \f$\dot{x}\f$ of an element.
+     * @param dxdot_dx        Value of \f$\dot{x} \cdot dx\f$.
+     * @param dx_dx           Value of \f$ x \cdot dx\f$.
+     * @param local_M_data    Mass matrix of an element, which takes the form of
+     *                        \f$ \inta N^T N\mathrm{d}\Omega\f$. Not used.
+     * @param local_K_data    Laplacian matrix of an element, which takes the
+     *         form of \f$ \int (\nabla N)^T K \nabla N\mathrm{d}\Omega\f$.
+     *                        Not used.
+     * @param local_b_data    Right hand side vector of an element.
+     * @param local_Jac_data  Element Jacobian matrix for the Newton-Raphson
+     *                        method.
+     * @param local_coupled_solutions Nodal values of solutions of the coupled
+     *                                processes of an element.
+     */
+    void assembleWithJacobianForPressureEquations(
+        double const t, std::vector<double> const& local_xdot,
+        const double dxdot_dx, const double dx_dx,
+        std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+        std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
+        LocalCoupledSolutions const& local_coupled_solutions);
+
+private:
+    RichardsMechanicsProcessData<DisplacementDim>& _process_data;
+
+    using BMatricesType =
+        BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+    using IpData =
+        IntegrationPointData<BMatricesType, ShapeMatricesTypeDisplacement,
+                             ShapeMatricesTypePressure, DisplacementDim,
+                             ShapeFunctionDisplacement::NPOINTS>;
+    std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+
+    IntegrationMethod _integration_method;
+    MeshLib::Element const& _element;
+    bool const _is_axially_symmetric;
+    SecondaryData<
+        typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
+        _secondary_data;
+
+    static const int pressure_index = 0;
+    static const int pressure_size = ShapeFunctionPressure::NPOINTS;
+    static const int displacement_index = ShapeFunctionPressure::NPOINTS;
+    static const int displacement_size =
+        ShapeFunctionDisplacement::NPOINTS * DisplacementDim;
+};
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib
+
+#include "RichardsMechanicsFEM-impl.h"

ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "RichardsMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MeshLib/Elements/Utils.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/RichardsMechanics/CreateLocalAssemblers.h"
+
+#include "RichardsMechanicsFEM.h"
+#include "RichardsMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <int DisplacementDim>
+RichardsMechanicsProcess<DisplacementDim>::RichardsMechanicsProcess(
+    MeshLib::Mesh& mesh,
+    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order,
+    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+        process_variables,
+    RichardsMechanicsProcessData<DisplacementDim>&& process_data,
+    SecondaryVariableCollection&& secondary_variables,
+    NumLib::NamedFunctionCaller&& named_function_caller,
+    bool const use_monolithic_scheme)
+    : Process(mesh, std::move(jacobian_assembler), parameters,
+              integration_order, std::move(process_variables),
+              std::move(secondary_variables), std::move(named_function_caller),
+              use_monolithic_scheme),
+      _process_data(std::move(process_data))
+{
+    _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+    _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
+        mesh, "HydraulicFlow", MeshLib::MeshItemType::Node, 1);
+}
+
+template <int DisplacementDim>
+bool RichardsMechanicsProcess<DisplacementDim>::isLinear() const
+{
+    return false;
+}
+
+template <int DisplacementDim>
+MathLib::MatrixSpecifications
+RichardsMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
+    const int process_id) const
+{
+    // For the monolithic scheme or the M process (deformation) in the staggered
+    // scheme.
+    if (_use_monolithic_scheme || process_id == 1)
+    {
+        auto const& l = *_local_to_global_index_map;
+        return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+                &l.getGhostIndices(), &this->_sparsity_pattern};
+    }
+
+    // For staggered scheme and H process (pressure).
+    auto const& l = *_local_to_global_index_map_with_base_nodes;
+    return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+            &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::constructDofTable()
+{
+    // Create single component dof in every of the mesh's nodes.
+    _mesh_subset_all_nodes =
+        std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
+    // Create single component dof in the mesh's base nodes.
+    _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
+    _mesh_subset_base_nodes =
+        std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
+
+    // TODO move the two data members somewhere else.
+    // for extrapolation of secondary variables of stress or strain
+    std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
+        *_mesh_subset_all_nodes};
+    _local_to_global_index_map_single_component =
+        std::make_unique<NumLib::LocalToGlobalIndexMap>(
+            std::move(all_mesh_subsets_single_component),
+            // by location order is needed for output
+            NumLib::ComponentOrder::BY_LOCATION);
+
+    if (_use_monolithic_scheme)
+    {
+        // For pressure, which is the first
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets{
+            *_mesh_subset_base_nodes};
+
+        // For displacement.
+        const int monolithic_process_id = 0;
+        std::generate_n(std::back_inserter(all_mesh_subsets),
+                        getProcessVariables(monolithic_process_id)[1]
+                            .get()
+                            .getNumberOfComponents(),
+                        [&]() { return *_mesh_subset_all_nodes; });
+
+        std::vector<int> const vec_n_components{1, DisplacementDim};
+        _local_to_global_index_map =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets), vec_n_components,
+                NumLib::ComponentOrder::BY_LOCATION);
+        assert(_local_to_global_index_map);
+    }
+    else
+    {
+        // For displacement equation.
+        const int process_id = 1;
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets;
+        std::generate_n(
+            std::back_inserter(all_mesh_subsets),
+            getProcessVariables(process_id)[0].get().getNumberOfComponents(),
+            [&]() { return *_mesh_subset_all_nodes; });
+
+        std::vector<int> const vec_n_components{DisplacementDim};
+        _local_to_global_index_map =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets), vec_n_components,
+                NumLib::ComponentOrder::BY_LOCATION);
+
+        // For pressure equation.
+        // Collect the mesh subsets with base nodes in a vector.
+        std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
+            *_mesh_subset_base_nodes};
+        _local_to_global_index_map_with_base_nodes =
+            std::make_unique<NumLib::LocalToGlobalIndexMap>(
+                std::move(all_mesh_subsets_base_nodes),
+                // by location order is needed for output
+                NumLib::ComponentOrder::BY_LOCATION);
+
+        _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
+            *_local_to_global_index_map_with_base_nodes, _mesh);
+
+        assert(_local_to_global_index_map);
+        assert(_local_to_global_index_map_with_base_nodes);
+    }
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    MeshLib::Mesh const& mesh,
+    unsigned const integration_order)
+{
+    const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
+    const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
+    ProcessLib::RichardsMechanics::createLocalAssemblers<
+        DisplacementDim, RichardsMechanicsLocalAssembler>(
+        mesh.getDimension(), mesh.getElements(), dof_table,
+        // use displacement process variable to set shape function order
+        getProcessVariables(mechanical_process_id)[deformation_variable_id]
+            .get()
+            .getShapeFunctionOrder(),
+        _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+        _process_data);
+
+    _secondary_variables.addSecondaryVariable(
+        "sigma",
+        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+                             DisplacementDim>::RowsAtCompileTime,
+                         getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtSigma));
+
+    _secondary_variables.addSecondaryVariable(
+        "epsilon",
+        makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
+                             DisplacementDim>::RowsAtCompileTime,
+                         getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtEpsilon));
+
+    _secondary_variables.addSecondaryVariable(
+        "velocity",
+        makeExtrapolator(DisplacementDim, getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtDarcyVelocity));
+
+    _secondary_variables.addSecondaryVariable(
+        "saturation",
+        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+                         &LocalAssemblerInterface::getIntPtSaturation));
+
+    _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
+        MeshLib::MeshItemType::Cell, 1);
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
+{
+    if (_use_monolithic_scheme)
+    {
+        const int monolithic_process_id = 0;
+        initializeProcessBoundaryConditionsAndSourceTerms(
+            *_local_to_global_index_map, monolithic_process_id);
+        return;
+    }
+
+    // Staggered scheme:
+    // for the equations of pressure
+    const int hydraulic_process_id = 0;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
+
+    // for the equations of deformation.
+    const int mechanical_process_id = 1;
+    initializeProcessBoundaryConditionsAndSourceTerms(
+        *_local_to_global_index_map, mechanical_process_id);
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
+    const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+    GlobalVector& b)
+{
+    DBUG("Assemble the equations for RichardsMechanics");
+
+    // Note: This assembly function is for the Picard nonlinear solver. Since
+    // only the Newton-Raphson method is employed to simulate coupled HM
+    // processes in this class, this function is actually not used so far.
+
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_table = {std::ref(*_local_to_global_index_map)};
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+        dof_table, t, x, M, K, b, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::
+    assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
+                                        GlobalVector const& xdot,
+                                        const double dxdot_dx,
+                                        const double dx_dx, GlobalMatrix& M,
+                                        GlobalMatrix& K, GlobalVector& b,
+                                        GlobalMatrix& Jac)
+{
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    // For the monolithic scheme
+    if (_use_monolithic_scheme)
+    {
+        DBUG(
+            "Assemble the Jacobian of RichardsMechanics for the monolithic"
+            " scheme.");
+        dof_tables.emplace_back(*_local_to_global_index_map);
+    }
+    else
+    {
+        // For the staggered scheme
+        if (_coupled_solutions->process_id == 0)
+        {
+            DBUG(
+                "Assemble the Jacobian equations of liquid fluid process in "
+                "RichardsMechanics for the staggered scheme.");
+        }
+        else
+        {
+            DBUG(
+                "Assemble the Jacobian equations of mechanical process in "
+                "RichardsMechanics for the staggered scheme.");
+        }
+        dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
+        dof_tables.emplace_back(*_local_to_global_index_map);
+    }
+
+    GlobalExecutor::executeMemberDereferenced(
+        _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+        _local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
+        Jac, _coupled_solutions);
+
+    auto copyRhs = [&](int const variable_id, auto& output_vector) {
+        if (_use_monolithic_scheme)
+        {
+            transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
+                                              output_vector,
+                                              std::negate<double>());
+        }
+        else
+        {
+            transformVariableFromGlobalVector(
+                b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
+                std::negate<double>());
+        }
+    };
+    if (_use_monolithic_scheme || _coupled_solutions->process_id == 0)
+    {
+        copyRhs(0, *_hydraulic_flow);
+    }
+    if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
+    {
+        copyRhs(1, *_nodal_forces);
+    }
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
+    GlobalVector const& x, double const t, double const dt,
+    const int process_id)
+{
+    DBUG("PreTimestep RichardsMechanicsProcess.");
+
+    _process_data.dt = dt;
+    _process_data.t = t;
+
+    if (hasMechanicalProcess(process_id))
+        GlobalExecutor::executeMemberOnDereferenced(
+            &LocalAssemblerInterface::preTimestep, _local_assemblers,
+            *_local_to_global_index_map, x, t, dt);
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<
+    DisplacementDim>::postNonLinearSolverConcreteProcess(GlobalVector const& x,
+                                                         const double t,
+                                                         const int process_id)
+{
+    if (!hasMechanicalProcess(process_id))
+    {
+        return;
+    }
+
+    DBUG("PostNonLinearSolver RichardsMechanicsProcess.");
+    // Calculate strain, stress or other internal variables of mechanics.
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
+        getDOFTable(process_id), x, t, _use_monolithic_scheme);
+}
+
+template <int DisplacementDim>
+void RichardsMechanicsProcess<
+    DisplacementDim>::computeSecondaryVariableConcrete(const double t,
+                                                       GlobalVector const& x)
+{
+    DBUG("Compute the secondary variables for RichardsMechanicsProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+        *_local_to_global_index_map, t, x, _coupled_solutions);
+}
+
+template <int DisplacementDim>
+std::tuple<NumLib::LocalToGlobalIndexMap*, bool> RichardsMechanicsProcess<
+    DisplacementDim>::getDOFTableForExtrapolatorData() const
+{
+    const bool manage_storage = false;
+    return std::make_tuple(_local_to_global_index_map_single_component.get(),
+                           manage_storage);
+}
+
+template <int DisplacementDim>
+NumLib::LocalToGlobalIndexMap const&
+RichardsMechanicsProcess<DisplacementDim>::getDOFTable(
+    const int process_id) const
+{
+    if (hasMechanicalProcess(process_id))
+    {
+        return *_local_to_global_index_map;
+    }
+
+    // For the equation of pressure
+    return *_local_to_global_index_map_with_base_nodes;
+}
+
+template class RichardsMechanicsProcess<2>;
+template class RichardsMechanicsProcess<3>;
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "LocalAssemblerInterface.h"
+#include "ProcessLib/Process.h"
+#include "RichardsMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+struct LocalAssemblerInterface;
+
+/// Linear kinematics poro-mechanical/biphasic (fluid-solid mixture) model.
+///
+/// The mixture momentum balance and the mixture mass balance are solved under
+/// fully saturated conditions.
+template <int DisplacementDim>
+class RichardsMechanicsProcess final : public Process
+{
+public:
+    RichardsMechanicsProcess(
+        MeshLib::Mesh& mesh,
+        std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+            jacobian_assembler,
+        std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+        unsigned const integration_order,
+        std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+            process_variables,
+        RichardsMechanicsProcessData<DisplacementDim>&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller,
+        bool const use_monolithic_scheme);
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override;
+    //! @}
+
+    /**
+     * Get the size and the sparse pattern of the global matrix in order to
+     * create the global matrices and vectors for the system equations of this
+     * process.
+     *
+     * @param process_id Process ID. If the monolithic scheme is applied,
+     *                               process_id = 0. For the staggered scheme,
+     *                               process_id = 0 represents the
+     *                               hydraulic (H) process, while process_id = 1
+     *                               represents the mechanical (M) process.
+     * @return Matrix specifications including size and sparse pattern.
+     */
+    MathLib::MatrixSpecifications getMatrixSpecifications(
+        const int process_id) const override;
+
+private:
+    void constructDofTable() override;
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override;
+
+    void initializeBoundaryConditions() override;
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override;
+
+    void assembleWithJacobianConcreteProcess(
+        const double t, GlobalVector const& x, GlobalVector const& xdot,
+        const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+        GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+                                    double const dt,
+                                    const int process_id) override;
+
+    void postNonLinearSolverConcreteProcess(GlobalVector const& x,
+                                            const double t,
+                                            int const process_id) override;
+
+    NumLib::LocalToGlobalIndexMap const& getDOFTable(
+        const int process_id) const override;
+
+private:
+    std::vector<MeshLib::Node*> _base_nodes;
+    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
+    RichardsMechanicsProcessData<DisplacementDim> _process_data;
+
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_single_component;
+
+    /// Local to global index mapping for base nodes, which is used for linear
+    /// interpolation for pressure in the staggered scheme.
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
+        _local_to_global_index_map_with_base_nodes;
+
+    /// Sparsity pattern for the flow equation, and it is initialized only if
+    /// the staggered scheme is used.
+    GlobalSparsityPattern _sparsity_pattern_with_linear_element;
+
+    /// Solutions of the previous time step
+    std::array<std::unique_ptr<GlobalVector>, 2> _xs_previous_timestep;
+
+    void computeSecondaryVariableConcrete(const double t,
+                                          GlobalVector const& x) override;
+    /**
+     * @copydoc ProcessLib::Process::getDOFTableForExtrapolatorData()
+     */
+    std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
+    getDOFTableForExtrapolatorData() const override;
+
+    /// Check whether the process represented by \c process_id is/has
+    /// mechanical process. In the present implementation, the mechanical
+    /// process has process_id == 1 in the staggered scheme.
+    bool hasMechanicalProcess(int const process_id) const
+    {
+        return _use_monolithic_scheme || process_id == 1;
+    }
+
+    MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+    MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
+};
+
+extern template class RichardsMechanicsProcess<2>;
+extern template class RichardsMechanicsProcess<3>;
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib

ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include <memory>
+#include <utility>
+
+#include <Eigen/Dense>
+
+#include "ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+}
+}  // namespace MaterialLib
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <int DisplacementDim>
+struct RichardsMechanicsProcessData
+{
+    RichardsMechanicsProcessData(
+        std::unique_ptr<
+            ProcessLib::RichardsFlow::RichardsFlowMaterialProperties>&&
+            flow_material_,
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
+            solid_material_,
+        Parameter<double> const& intrinsic_permeability_,
+        Parameter<double> const& fluid_bulk_modulus_,
+        Parameter<double> const& biot_coefficient_,
+        Parameter<double> const& solid_density_,
+        Parameter<double> const& solid_bulk_modulus_,
+        Parameter<double> const& temperature_,
+        Eigen::Matrix<double, DisplacementDim, 1>
+            specific_body_force_)
+        : flow_material{std::move(flow_material_)},
+          solid_material{std::move(solid_material_)},
+          intrinsic_permeability(intrinsic_permeability_),
+          fluid_bulk_modulus(fluid_bulk_modulus_),
+          biot_coefficient(biot_coefficient_),
+          solid_density(solid_density_),
+          solid_bulk_modulus(solid_bulk_modulus_),
+          temperature(temperature_),
+          specific_body_force(std::move(specific_body_force_))
+    {
+    }
+
+    RichardsMechanicsProcessData(RichardsMechanicsProcessData&& other) =
+        default;
+
+    //! Copies are forbidden.
+    RichardsMechanicsProcessData(RichardsMechanicsProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(RichardsMechanicsProcessData const&) = delete;
+    void operator=(RichardsMechanicsProcessData&&) = delete;
+
+    std::unique_ptr<ProcessLib::RichardsFlow::RichardsFlowMaterialProperties>
+        flow_material;
+
+    /// The constitutive relation for the mechanical part.
+    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+        solid_material;
+    /// Permeability of the solid. A scalar quantity, Parameter<double>.
+    Parameter<double> const& intrinsic_permeability;
+    /// Fluid's bulk modulus. A scalar quantity, Parameter<double>.
+    Parameter<double> const& fluid_bulk_modulus;
+    /// Biot coefficient. A scalar quantity, Parameter<double>.
+    Parameter<double> const& biot_coefficient;
+    /// Density of the solid. A scalar quantity, Parameter<double>.
+    Parameter<double> const& solid_density;
+    /// Solid's bulk modulus. A scalar quantity, Parameter<double>.
+    Parameter<double> const& solid_bulk_modulus;
+    /// Reference temperature for material properties. A scalar quantity,
+    /// Parameter<double>.
+    Parameter<double> const& temperature;
+    /// Specific body forces applied to solid and fluid.
+    /// It is usually used to apply gravitational forces.
+    /// A vector of displacement dimension's length.
+    Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+    double dt = 0.0;
+    double t = 0.0;
+
+    MeshLib::PropertyVector<double>* element_saturation = nullptr;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+}  // namespace RichardsMechanics
+}  // namespace ProcessLib