[PL] Implemention of HydroMechanics process.

ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "CreateHydroMechanicsProcess.h"
+
+#include <cassert>
+
+#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+
+#include "HydroMechanicsProcess.h"
+#include "HydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+template <int DisplacementDim>
+class HydroMechanicsProcess;
+
+extern template class HydroMechanicsProcess<2>;
+extern template class HydroMechanicsProcess<3>;
+
+template <int DisplacementDim>
+std::unique_ptr<Process> createHydroMechanicsProcess(
+    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{process__type}
+    config.checkConfigParameter("type", "HYDRO_MECHANICS");
+    DBUG("Create HydroMechanicsProcess.");
+
+    // Process variable.
+    auto process_variables = findProcessVariables(
+        variables, config,
+        {//! \ogs_file_param_special{process__HYDRO_MECHANICS_process_variables__process_variable}
+          "pressure", "displacement"});
+
+    DBUG("Associate displacement with process variable \'%s\'.",
+         process_variables[1].get().getName().c_str());
+
+    if (process_variables[1].get().getNumberOfComponents() !=
+        DisplacementDim)
+    {
+        OGS_FATAL(
+            "Number of components of the process variable '%s' is different "
+            "from the displacement dimension: got %d, expected %d",
+            process_variables[1].get().getName().c_str(),
+            process_variables[1].get().getNumberOfComponents(),
+            DisplacementDim);
+    }
+
+    DBUG("Associate pressure with process variable \'%s\'.",
+         process_variables[0].get().getName().c_str());
+    if (process_variables[0].get().getNumberOfComponents() != 1)
+    {
+        OGS_FATAL(
+            "Pressure process variable '%s' is not a scalar variable but has "
+            "%d components.",
+            process_variables[0].get().getName().c_str(),
+            process_variables[0].get().getNumberOfComponents());
+    }
+
+
+    // Constitutive relation.
+    // read type;
+    auto const constitutive_relation_config =
+        //! \ogs_file_param{process__HYDRO_MECHANICS_constitutive_relation}
+        config.getConfigSubtree("constitutive_relation");
+
+    auto const type =
+        constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+        material = nullptr;
+    if (type == "LinearElasticIsotropic")
+    {
+        material =
+            MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
+                parameters, constitutive_relation_config);
+    }
+    else
+    {
+        OGS_FATAL(
+            "Cannot construct constitutive relation of given type \'%s\'.",
+            type.c_str());
+    }
+
+    // Intrinsic permeability
+    auto& intrinsic_permeability = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_intrinsic_permeability}
+        "intrinsic_permeability",
+        parameters, 1);
+
+    DBUG("Use \'%s\' as intrinsic conductivity parameter.",
+         intrinsic_permeability.name.c_str());
+
+    // Storage coefficient
+    auto& specific_storage = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_specific_storage}
+        "specific_storage", parameters, 1);
+
+    DBUG("Use \'%s\' as storage coefficient parameter.",
+         specific_storage.name.c_str());
+
+    // Fluid viscosity
+    auto& fluid_viscosity = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_fluid_viscosity}
+        "fluid_viscosity",
+        parameters, 1);
+    DBUG("Use \'%s\' as fluid viscosity parameter.",
+         fluid_viscosity.name.c_str());
+
+    // Fluid density
+    auto& fluid_density = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_fluid_density}
+        "fluid_density",
+        parameters, 1);
+    DBUG("Use \'%s\' as fluid density parameter.",
+         fluid_density.name.c_str());
+
+    // Biot coefficient
+    auto& biot_coefficient = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_biot_coefficient}
+        "biot_coefficient",
+        parameters, 1);
+    DBUG("Use \'%s\' as Biot coefficient parameter.",
+         biot_coefficient.name.c_str());
+
+    // Porosity
+    auto& porosity = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_porosity}
+        "porosity",
+        parameters, 1);
+    DBUG("Use \'%s\' as porosity parameter.",
+         porosity.name.c_str());
+
+    // Solid density
+    auto& solid_density = findParameter<double>(
+        config,
+        //! \ogs_file_param_special{process__HYDRO_MECHANICS_solid_density}
+        "solid_density",
+        parameters, 1);
+    DBUG("Use \'%s\' as solid density parameter.",
+         solid_density.name.c_str());
+
+    // Specific body force
+    Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+    {
+        std::vector<double> const b =
+            //! \ogs_file_param_special{process__HYDRO_MECHANICS_specific_body_force}
+            config.getConfigParameter<std::vector<double>>(
+                "specific_body_force");
+        if (specific_body_force.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",
+                specific_body_force.size(), DisplacementDim);
+
+        std::copy_n(b.data(), b.size(), specific_body_force.data());
+    }
+
+    HydroMechanicsProcessData<DisplacementDim> process_data{
+        std::move(material),
+        intrinsic_permeability,
+        specific_storage,
+        fluid_viscosity,
+        fluid_density,
+        biot_coefficient,
+        porosity,
+        solid_density,
+        specific_body_force};
+
+    SecondaryVariableCollection secondary_variables;
+
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"HydroMechanics_displacement"});
+
+    ProcessLib::parseSecondaryVariables(config, secondary_variables,
+                                        named_function_caller);
+
+    return std::unique_ptr<HydroMechanicsProcess<DisplacementDim>>{
+        new HydroMechanicsProcess<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)}};
+}
+
+template std::unique_ptr<Process> createHydroMechanicsProcess<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> createHydroMechanicsProcess<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 HydroMechanics
+}  // namespace ProcessLib

ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_CREATEHYDROMECHANICSPROCESS_H_
+#define PROCESS_LIB_CREATEHYDROMECHANICSPROCESS_H_
+
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createHydroMechanicsProcess(
+    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 HydroMechanics
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_CREATEHYDROMECHANICSPROCESS_H_

ProcessLib/HydroMechanics/CreateLocalAssemblers.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_HYDROMECHANICS_CREATE_LOCAL_ASSEMBLERS_H_
+#define PROCESSLIB_HYDROMECHANICS_CREATE_LOCAL_ASSEMBLERS_H_
+
+#include <vector>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+#include "LocalDataInitializer.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+namespace detail
+{
+template <unsigned GlobalDim, int DisplacementDim,
+          template <typename, typename, typename, unsigned, 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,
+                             DisplacementDim, 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 DisplacementDim, template <typename, typename, typename, unsigned, 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.");
+
+    switch (dimension)
+    {
+        case 2:
+            detail::createLocalAssemblers<2, DisplacementDim,
+                                          LocalAssemblerImplementation>(
+                dof_table, shapefunction_order, mesh_elements, local_assemblers,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
+        case 3:
+            detail::createLocalAssemblers<3, DisplacementDim,
+                                          LocalAssemblerImplementation>(
+                dof_table, shapefunction_order, mesh_elements, local_assemblers,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
+        default:
+            OGS_FATAL(
+                "Meshes with dimension different than two and three are not "
+                "supported.");
+    }
+}
+}  // HydroMechanics
+
+}  // ProcessLib
+
+#endif  // PROCESSLIB_HYDROMECHANICS_CREATE_LOCAL_ASSEMBLERS_H_

ProcessLib/HydroMechanics/HydroMechanicsFEM.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HYDROMECHANICS_FEM_H_
+#define PROCESS_LIB_HYDROMECHANICS_FEM_H_
+
+#include <iostream>
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/SolidModels/KelvinVector.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "HydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+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())
+    {
+    }
+
+#if defined(_MSC_VER) && _MSC_VER < 1900
+    // The default generated move-ctor is correctly generated for other
+    // compilers.
+    explicit IntegrationPointData(IntegrationPointData&& other)
+        : b_matrices(std::move(other.b_matrices)),
+          sigma_eff(std::move(other.sigma_eff)),
+          sigma_eff_prev(std::move(other.sigma_eff_prev)),
+          eps(std::move(other.eps)),
+          eps_prev(std::move(other.eps_prev)),
+          solid_material(other.solid_material),
+          material_state_variables(std::move(other.material_state_variables)),
+          C(std::move(other.C)),
+          integration_weight(std::move(other.integration_weight))
+    {
+    }
+#endif  // _MSC_VER
+
+    typename ShapeMatrixTypeDisplacement::template MatrixType<
+        DisplacementDim, NPoints * DisplacementDim>
+        N_u;
+    //typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
+    typename BMatricesType::BMatrixType b_matrices;
+    typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
+    typename BMatricesType::KelvinVectorType eps, eps_prev;
+
+    typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
+    typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
+
+    MaterialLib::Solids::MechanicsBase<DisplacementDim>& solid_material;
+    std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
+        DisplacementDim>::MaterialStateVariables>
+        material_state_variables;
+
+    typename BMatricesType::KelvinMatrixType C;
+    double integration_weight;
+
+    void pushBackState()
+    {
+        eps_prev = eps;
+        sigma_eff_prev = sigma_eff;
+        material_state_variables->pushBackState();
+    }
+
+    template <typename DisplacementVectorType>
+    void updateConstitutiveRelation(double const t,
+                                    SpatialPosition const& x_position,
+                                    double const dt,
+                                    DisplacementVectorType const& u)
+    {
+        eps.noalias() = b_matrices * u;
+        solid_material.computeConstitutiveRelation(
+            t, x_position, dt, eps_prev, eps, sigma_eff_prev, sigma_eff, C,
+            *material_state_variables);
+    }
+};
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int DisplacementDim>
+class HydroMechanicsLocalAssembler : public ProcessLib::LocalAssemblerInterface
+{
+public:
+    using ShapeMatricesTypeDisplacement =
+        ShapeMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+
+    // Types for pressure.
+    using ShapeMatricesTypePressure =
+        ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
+
+    HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete;
+    HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete;
+
+    HydroMechanicsLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const /*local_matrix_size*/,
+        bool is_axially_symmetric,
+        unsigned const integration_order,
+        HydroMechanicsProcessData<DisplacementDim>& process_data)
+        : _process_data(process_data),
+          _integration_method(integration_order),
+          _element(e)
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        _ip_data.reserve(n_integration_points);
+
+        auto const shape_matrices_u =
+            initShapeMatrices<ShapeFunctionDisplacement,
+                              ShapeMatricesTypeDisplacement, IntegrationMethod,
+                              DisplacementDim>(e, is_axially_symmetric,
+                                               _integration_method);
+
+        auto const shape_matrices_p =
+            initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
+                              IntegrationMethod, DisplacementDim>(
+                e, is_axially_symmetric, _integration_method);
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            // displacement (subscript u)
+            _ip_data.emplace_back(*_process_data.material);
+            auto& ip_data = _ip_data[ip];
+            _ip_data[ip].integration_weight =
+                _integration_method.getWeightedPoint(ip).getWeight() *
+                shape_matrices_u[ip].detJ;
+            ip_data.b_matrices.resize(
+                kelvin_vector_size,
+                ShapeFunctionDisplacement::NPOINTS * DisplacementDim);
+
+            auto const x_coord =
+                interpolateXCoordinate<ShapeFunctionDisplacement,
+                                       ShapeMatricesTypeDisplacement>(
+                    e, shape_matrices_u[ip].N);
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS>(
+                shape_matrices_u[ip].dNdx, ip_data.b_matrices,
+                is_axially_symmetric, shape_matrices_u[ip].N, x_coord);
+
+            ip_data.sigma_eff.resize(kelvin_vector_size);
+            ip_data.sigma_eff_prev.resize(kelvin_vector_size);
+            ip_data.eps.resize(kelvin_vector_size);
+            ip_data.eps_prev.resize(kelvin_vector_size);
+            ip_data.C.resize(kelvin_vector_size, kelvin_vector_size);
+
+            ip_data.N_u = ShapeMatricesTypeDisplacement::template MatrixType<
+                DisplacementDim, displacement_size>::Zero(DisplacementDim,
+                                                          displacement_size);
+            for (int i = 0; i < DisplacementDim; ++i)
+                ip_data.N_u
+                    .template block<1, displacement_size / DisplacementDim>(
+                        i, i * displacement_size / DisplacementDim)
+                    .noalias() = shape_matrices_u[ip].N;
+
+            ip_data.N_p = shape_matrices_p[ip].N;
+            ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
+        }
+    }
+
+    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
+    {
+        OGS_FATAL(
+            "HydroMechanicsLocalAssembler: assembly without jacobian is not "
+            "implemented.");
+    }
+
+    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
+    {
+        assert(local_x.size() == pressure_size + displacement_size);
+
+        auto p =
+            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+                pressure_size> const>(local_x.data() + pressure_index,
+                                      pressure_size);
+
+        auto u = Eigen::Map<typename ShapeMatricesTypeDisplacement::
+                                template VectorType<displacement_size> const>(
+            local_x.data() + displacement_index, displacement_size);
+
+        auto p_dot =
+            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+                pressure_size> const>(local_xdot.data() + pressure_index,
+                                      pressure_size);
+        auto u_dot =
+            Eigen::Map<typename ShapeMatricesTypeDisplacement::
+                           template VectorType<displacement_size> const>(
+                local_xdot.data() + displacement_index, displacement_size);
+
+        auto local_Jac = MathLib::createZeroedMatrix<
+            typename ShapeMatricesTypeDisplacement::template MatrixType<
+                displacement_size + pressure_size,
+                displacement_size + pressure_size>>(
+            local_Jac_data, displacement_size + pressure_size,
+            displacement_size + pressure_size);
+
+        auto local_rhs = MathLib::createZeroedVector<
+            typename ShapeMatricesTypeDisplacement::template VectorType<
+                displacement_size + pressure_size>>(
+            local_rhs_data, displacement_size + pressure_size);
+
+        typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
+        laplace_p.setZero(pressure_size, pressure_size);
+
+        typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
+        storage_p.setZero(pressure_size, pressure_size);
+
+        typename ShapeMatricesTypeDisplacement::template MatrixType<
+            displacement_size, pressure_size>
+            Kup;
+        Kup.setZero(displacement_size, pressure_size);
+
+        double const& dt = _process_data.dt;
+
+        SpatialPosition x_position;
+        x_position.setElementID(_element.getID());
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            x_position.setIntegrationPoint(ip);
+            auto const& w = _ip_data[ip].integration_weight;
+
+            auto const& N_p = _ip_data[ip].N_p;
+            auto const& N_u = _ip_data[ip].N_u;
+            auto const& dNdx_p = _ip_data[ip].dNdx_p;
+
+            auto const& B = _ip_data[ip].b_matrices;
+            auto const& sigma_eff = _ip_data[ip].sigma_eff;
+
+            auto const& C = _ip_data[ip].C;
+
+            double const S =
+                _process_data.specific_storage(t, x_position)[0];
+            double const K_over_mu =
+                _process_data.intrinsic_permeability(t, x_position)[0] /
+                _process_data.fluid_viscosity(t, x_position)[0];
+            auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
+            auto const rho_sr = _process_data.solid_density(t, x_position)[0];
+            auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+            auto const porosity = _process_data.porosity(t, x_position)[0];
+            auto const& b = _process_data.specific_body_force;
+            auto const& identity2 = MaterialLib::SolidModels::Invariants<
+                kelvin_vector_size>::identity2;
+
+            //
+            // displacement equation, displacement part
+            //
+            _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u);
+
+            local_Jac
+                .template block<displacement_size, displacement_size>(
+                    displacement_index, displacement_index)
+                .noalias() += B.transpose() * C * B * w;
+
+            double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
+            local_rhs.template segment<displacement_size>(displacement_index)
+                .noalias() -=
+                (B.transpose() * sigma_eff - N_u.transpose() * rho * b) * w;
+
+            //
+            // displacement equation, pressure part
+            //
+            Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
+
+            //
+            // pressure equation, pressure part.
+            //
+            laplace_p.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
+
+            storage_p.noalias() += N_p.transpose() * S * N_p * w;
+
+            local_rhs.template segment<pressure_size>(pressure_index)
+                .noalias() += dNdx_p.transpose() * rho_fr * K_over_mu * b * w;
+
+            //
+            // pressure equation, displacement part.
+            //
+            // Reusing Kup.transpose().
+        }
+        // displacement equation, pressure part
+        local_Jac
+            .template block<displacement_size, pressure_size>(
+                displacement_index, pressure_index)
+            .noalias() -= Kup;
+
+        // pressure equation, pressure part.
+        local_Jac
+            .template block<pressure_size, pressure_size>(pressure_index,
+                                                          pressure_index)
+            .noalias() += laplace_p + storage_p / dt;
+
+        // pressure equation, displacement part.
+        local_Jac
+            .template block<pressure_size, displacement_size>(
+                pressure_index, displacement_index)
+            .noalias() += Kup.transpose() / dt;
+
+        // pressure equation
+        local_rhs.template segment<pressure_size>(pressure_index)
+            .noalias() -=
+            laplace_p * p + storage_p * p_dot + Kup.transpose() * u_dot;
+
+        // displacement equation
+        local_rhs.template segment<displacement_size>(displacement_index)
+            .noalias() += Kup * p;
+    }
+
+    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();
+        }
+    }
+
+private:
+    HydroMechanicsProcessData<DisplacementDim>& _process_data;
+
+    using BMatricesType =
+        BMatrixPolicyType<ShapeFunctionDisplacement, DisplacementDim>;
+    std::vector<IntegrationPointData<
+        BMatricesType, ShapeMatricesTypeDisplacement, ShapeMatricesTypePressure,
+        DisplacementDim, ShapeFunctionDisplacement::NPOINTS>>
+        _ip_data;
+
+    IntegrationMethod _integration_method;
+    MeshLib::Element const& _element;
+
+    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;
+    static const int kelvin_vector_size =
+        KelvinVectorDimensions<DisplacementDim>::value;
+};
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+          typename IntegrationMethod, unsigned GlobalDim, int DisplacementDim>
+class LocalAssemblerData final
+    : public HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                          ShapeFunctionPressure,
+                                          IntegrationMethod, DisplacementDim>
+{
+public:
+    LocalAssemblerData(LocalAssemblerData const&) = delete;
+    LocalAssemblerData(LocalAssemblerData&&) = delete;
+
+    LocalAssemblerData(MeshLib::Element const& e,
+                       std::size_t const local_matrix_size,
+                       bool is_axially_symmetric,
+                       unsigned const integration_order,
+                       HydroMechanicsProcessData<DisplacementDim>& process_data)
+        : HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
+                                       ShapeFunctionPressure, IntegrationMethod,
+                                       DisplacementDim>(
+              e, local_matrix_size, is_axially_symmetric, integration_order,
+              process_data)
+    {
+    }
+};
+
+}  // namespace HydroMechanics
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_HYDROMECHANICS_FEM_H_

ProcessLib/HydroMechanics/HydroMechanicsProcess-fwd.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HYDROMECHANICSPROCESS_FWD_H_
+#define PROCESS_LIB_HYDROMECHANICSPROCESS_FWD_H_
+
+#include "HydroMechanicsProcess.h"
+
+extern template class ProcessLib::HydroMechanics::HydroMechanicsProcess<2>;
+extern template class ProcessLib::HydroMechanics::HydroMechanicsProcess<3>;
+
+#endif  // PROCESS_LIB_HYDROMECHANICSPROCESS_FWD_H_

ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "HydroMechanicsProcess-fwd.h"
+#include "HydroMechanicsProcess.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+
+template class HydroMechanicsProcess<2>;
+template class HydroMechanicsProcess<3>;
+
+}   // namespace HydroMechanics
+}   // namespace ProcessLib

ProcessLib/HydroMechanics/HydroMechanicsProcess.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_HYDROMECHANICSPROCESS_H_
+#define PROCESS_LIB_HYDROMECHANICSPROCESS_H_
+
+#include <cassert>
+
+#include "MeshLib/Elements/Utils.h"
+#include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
+#include "ProcessLib/Process.h"
+
+#include "HydroMechanicsFEM.h"
+#include "HydroMechanicsProcessData.h"
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+template <int DisplacementDim>
+class HydroMechanicsProcess final : public Process
+{
+    using Base = Process;
+
+public:
+    HydroMechanicsProcess(
+        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::reference_wrapper<ProcessVariable>>&&
+            process_variables,
+        HydroMechanicsProcessData<DisplacementDim>&& process_data,
+        SecondaryVariableCollection&& secondary_variables,
+        NumLib::NamedFunctionCaller&& named_function_caller)
+        : Process(mesh, std::move(jacobian_assembler), parameters,
+                  integration_order, std::move(process_variables),
+                  std::move(secondary_variables),
+                  std::move(named_function_caller)),
+          _process_data(std::move(process_data))
+    {
+    }
+
+    //! \name ODESystem interface
+    //! @{
+
+    bool isLinear() const override { return false; }
+    //! @}
+
+private:
+    void constructDofTable() override
+    {
+        // Create single component dof in every of the mesh's nodes.
+        _mesh_subset_all_nodes.reset(
+            new 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.reset(
+            new MeshLib::MeshSubset(_mesh, &_base_nodes));
+
+        // Collect the mesh subsets in a vector.
+
+        // For pressure, which is the first
+        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+        all_mesh_subsets.push_back(std::unique_ptr<MeshLib::MeshSubsets>{
+            new MeshLib::MeshSubsets{_mesh_subset_base_nodes.get()}});
+
+        // For displacement.
+        std::generate_n(
+            std::back_inserter(all_mesh_subsets),
+            getProcessVariables()[1].get().getNumberOfComponents(),
+            [&]() {
+                return std::unique_ptr<MeshLib::MeshSubsets>{
+                    new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+            });
+
+        std::vector<unsigned> const vec_n_components{1, DisplacementDim};
+        _local_to_global_index_map.reset(new NumLib::LocalToGlobalIndexMap(
+            std::move(all_mesh_subsets), vec_n_components,
+            NumLib::ComponentOrder::BY_LOCATION));
+    }
+
+    void initializeConcreteProcess(
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh,
+        unsigned const integration_order) override
+    {
+        ProcessLib::HydroMechanics::createLocalAssemblers<DisplacementDim,
+                                                          LocalAssemblerData>(
+            mesh.getDimension(), mesh.getElements(), dof_table,
+            // use displacment process variable for shapefunction order
+            getProcessVariables()[1].get().getShapeFunctionOrder(),
+            _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+            _process_data);
+    }
+
+    void assembleConcreteProcess(const double t, GlobalVector const& x,
+                                 GlobalMatrix& M, GlobalMatrix& K,
+                                 GlobalVector& b) override
+    {
+        DBUG("Assemble HydroMechanicsProcess.");
+
+        // Call global assembler for each local assembly item.
+        GlobalExecutor::executeMemberDereferenced(
+            _global_assembler, &VectorMatrixAssembler::assemble,
+            _local_assemblers, *_local_to_global_index_map, t, x, M, K, b);
+    }
+
+    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
+    {
+        DBUG("AssembleJacobian HydroMechanicsProcess.");
+
+        // Call global assembler for each local assembly item.
+        GlobalExecutor::executeMemberDereferenced(
+            _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+            _local_assemblers, *_local_to_global_index_map, t, x, xdot,
+            dxdot_dx, dx_dx, M, K, b, Jac);
+    }
+
+    void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+                                    double const dt) override
+    {
+        DBUG("PreTimestep HydroMechanicsProcess.");
+
+        _process_data.dt = dt;
+        _process_data.t = t;
+
+        GlobalExecutor::executeMemberOnDereferenced(
+            &LocalAssemblerInterface::preTimestep, _local_assemblers,
+            *_local_to_global_index_map, x, t, dt);
+    }
+
+    void postTimestepConcreteProcess(GlobalVector const& x) override
+    {
+        DBUG("PostTimestep HydroMechanicsProcess.");
+
+        GlobalExecutor::executeMemberOnDereferenced(
+            &LocalAssemblerInterface::postTimestep, _local_assemblers,
+            *_local_to_global_index_map, x);
+    }
+
+private:
+    std::vector<MeshLib::Node*> _base_nodes;
+    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
+    HydroMechanicsProcessData<DisplacementDim> _process_data;
+
+    std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+};
+
+}  // namespace HydroMechanics
+}  // namespace ProcessLib
+
+#endif  // PROCESS_LIB_HYDROMECHANICSPROCESS_H_

ProcessLib/HydroMechanics/HydroMechanicsProcessData.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_HYDROMECHANICS_HYDROMECHANICSPROCESSDATA_H_
+#define PROCESSLIB_HYDROMECHANICS_HYDROMECHANICSPROCESSDATA_H_
+
+#include <Eigen/Dense>
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+namespace HydroMechanics
+{
+template <int DisplacementDim>
+struct HydroMechanicsProcessData
+{
+    HydroMechanicsProcessData(
+        std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>&&
+            material_,
+        Parameter<double> const& intrinsic_permeability_,
+        Parameter<double> const& specific_storage_,
+        Parameter<double> const& fluid_viscosity_,
+        Parameter<double> const& fluid_density_,
+        Parameter<double> const& biot_coefficient_,
+        Parameter<double> const& porosity_,
+        Parameter<double> const& solid_density_,
+        Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
+        : material{std::move(material_)},
+          intrinsic_permeability(intrinsic_permeability_),
+          specific_storage(specific_storage_),
+          fluid_viscosity(fluid_viscosity_),
+          fluid_density(fluid_density_),
+          biot_coefficient(biot_coefficient_),
+          porosity(porosity_),
+          solid_density(solid_density_),
+          specific_body_force(specific_body_force_)
+    {
+    }
+
+    HydroMechanicsProcessData(HydroMechanicsProcessData&& other)
+        : material{std::move(other.material)},
+          intrinsic_permeability(other.intrinsic_permeability),
+          specific_storage(other.specific_storage),
+          fluid_viscosity(other.fluid_viscosity),
+          fluid_density(other.fluid_density),
+          biot_coefficient(other.biot_coefficient),
+          porosity(other.porosity),
+          solid_density(other.solid_density),
+          specific_body_force(other.specific_body_force),
+          dt(other.dt),
+          t(other.t)
+    {
+    }
+
+    //! Copies are forbidden.
+    HydroMechanicsProcessData(HydroMechanicsProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(HydroMechanicsProcessData const&) = delete;
+
+    //! Assignments are not needed.
+    void operator=(HydroMechanicsProcessData&&) = delete;
+
+    std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+        material;
+    Parameter<double> const& intrinsic_permeability;
+    Parameter<double> const& specific_storage;
+    Parameter<double> const& fluid_viscosity;
+    Parameter<double> const& fluid_density;
+    Parameter<double> const& biot_coefficient;
+    Parameter<double> const& porosity;
+    Parameter<double> const& solid_density;
+    Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+    double dt;
+    double t;
+};
+
+}  // namespace HydroMechanics
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_HYDROMECHANICS_HYDROMECHANICSPROCESSDATA_H_

ProcessLib/HydroMechanics/LocalDataInitializer.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_HYDROMECHANICS_LOCALDATAINITIALIZER_H_
+#define PROCESSLIB_HYDROMECHANICS_LOCALDATAINITIALIZER_H_
+
+#include <functional>
+#include <memory>
+#include <typeindex>
+#include <typeinfo>
+#include <type_traits>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/FiniteElement/LowerDimShapeTable.h"
+#include "NumLib/Fem/Integration/GaussIntegrationPolicy.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 HydroMechanics
+{
+/// 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, and provides additional template argument DisplacementDim.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, unsigned, int>
+          class LocalAssemblerData,
+          unsigned GlobalDim, int DisplacementDim, 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.",
+                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::GaussIntegrationPolicy<
+        typename ShapeFunctionDisplacement::MeshElement>::IntegrationMethod;
+
+    template <typename ShapeFunctionDisplacement,
+              typename ShapeFunctionPressure>
+    using LAData =
+        LocalAssemblerData<ShapeFunctionDisplacement, ShapeFunctionPressure,
+                           IntegrationMethod<ShapeFunctionDisplacement>,
+                           GlobalDim, DisplacementDim>;
+
+    /// 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 HydroMechanics
+}  // 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
+
+#endif  // PROCESSLIB_HYDROMECHANICS_LOCALDATAINITIALIZER_H_