From 2fa3356c73af3a7c43aa751a18935d069acfb745 Mon Sep 17 00:00:00 2001
From: Dmitrij Naumov <dmitrij@naumov.de>
Date: Wed, 27 Jun 2018 17:45:09 +0200
Subject: [PATCH] [PL] Implementation of RichardsMechanics process.
---
Applications/ApplicationsLib/ProjectData.cpp | 27 +
CMakeLists.txt | 1 +
ProcessLib/RichardsMechanics/CMakeLists.txt | 6 +
.../RichardsMechanics/CreateLocalAssemblers.h | 87 ++
.../CreateRichardsMechanicsProcess.cpp | 239 +++++
.../CreateRichardsMechanicsProcess.h | 65 ++
.../RichardsMechanics/IntegrationPointData.h | 96 ++
.../LocalAssemblerInterface.h | 48 +
.../RichardsMechanics/LocalDataInitializer.h | 289 ++++++
.../RichardsMechanicsFEM-impl.h | 833 ++++++++++++++++++
.../RichardsMechanics/RichardsMechanicsFEM.h | 237 +++++
.../RichardsMechanicsProcess.cpp | 381 ++++++++
.../RichardsMechanicsProcess.h | 140 +++
.../RichardsMechanicsProcessData.h | 104 +++
14 files changed, 2553 insertions(+)
create mode 100644 ProcessLib/RichardsMechanics/CMakeLists.txt
create mode 100644 ProcessLib/RichardsMechanics/CreateLocalAssemblers.h
create mode 100644 ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
create mode 100644 ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h
create mode 100644 ProcessLib/RichardsMechanics/IntegrationPointData.h
create mode 100644 ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
create mode 100644 ProcessLib/RichardsMechanics/LocalDataInitializer.h
create mode 100644 ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
create mode 100644 ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
create mode 100644 ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
create mode 100644 ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
create mode 100644 ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index d3d07d07bd6..10b42a97ca4 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/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")
{
diff --git a/CMakeLists.txt b/CMakeLists.txt
index add4d4668f4..ddd5670df00 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -99,6 +99,7 @@ set(ProcessesList
PhaseField
RichardsComponentTransport
RichardsFlow
+ RichardsMechanics
SmallDeformation
TES
ThermalTwoPhaseFlowWithPP
diff --git a/ProcessLib/RichardsMechanics/CMakeLists.txt b/ProcessLib/RichardsMechanics/CMakeLists.txt
new file mode 100644
index 00000000000..6e124260816
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/CMakeLists.txt
@@ -0,0 +1,6 @@
+APPEND_SOURCE_FILES(SOURCES)
+
+add_library(RichardsMechanics ${SOURCES})
+target_link_libraries(RichardsMechanics PUBLIC ProcessLib RichardsFlow)
+
+include(Tests.cmake)
diff --git a/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h b/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h
new file mode 100644
index 00000000000..7a7c825cfc0
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h
@@ -0,0 +1,87 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
new file mode 100644
index 00000000000..d901fb3f40a
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -0,0 +1,239 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h
new file mode 100644
index 00000000000..45c28c1b9a5
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.h
@@ -0,0 +1,65 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/IntegrationPointData.h b/ProcessLib/RichardsMechanics/IntegrationPointData.h
new file mode 100644
index 00000000000..cc03455f711
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/IntegrationPointData.h
@@ -0,0 +1,96 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
new file mode 100644
index 00000000000..d16c6081a3e
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
@@ -0,0 +1,48 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/LocalDataInitializer.h b/ProcessLib/RichardsMechanics/LocalDataInitializer.h
new file mode 100644
index 00000000000..f7cfc1dee1c
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/LocalDataInitializer.h
@@ -0,0 +1,289 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
new file mode 100644
index 00000000000..5d6f227360b
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -0,0 +1,833 @@
+/**
+ * \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
+ *
+ * \file RichardsMechanicsFEM-impl.h
+ * Created on November 29, 2017, 2:03 PM
+ */
+
+#pragma once
+
+#include "RichardsMechanicsFEM.h"
+
+#include "MathLib/KelvinVector.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
+
+namespace ProcessLib
+{
+namespace RichardsMechanics
+{
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ 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)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N_u.resize(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.solid_material);
+ auto& ip_data = _ip_data[ip];
+ auto const& sm_u = shape_matrices_u[ip];
+ _ip_data[ip].integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm_u.integralMeasure * sm_u.detJ;
+
+ ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
+ DisplacementDim, displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+ for (int i = 0; i < DisplacementDim; ++i)
+ ip_data.N_u_op
+ .template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = sm_u.N;
+
+ ip_data.N_u = sm_u.N;
+ ip_data.dNdx_u = sm_u.dNdx;
+
+ ip_data.N_p = shape_matrices_p[ip].N;
+ ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
+
+ _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::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)
+{
+ assert(local_x.size() == pressure_size + displacement_size);
+
+ auto p_L =
+ 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 K = MathLib::createZeroedMatrix<
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size + pressure_size,
+ displacement_size + pressure_size>>(
+ local_K_data, displacement_size + pressure_size,
+ displacement_size + pressure_size);
+
+ auto M = MathLib::createZeroedMatrix<
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size + pressure_size,
+ displacement_size + pressure_size>>(
+ local_M_data, displacement_size + pressure_size,
+ displacement_size + pressure_size);
+
+ auto rhs = MathLib::createZeroedVector<
+ typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size + pressure_size>>(
+ local_rhs_data, displacement_size + pressure_size);
+
+ double const& dt = _process_data.dt;
+ auto const material_id =
+ _process_data.flow_material->getMaterialID(_element.getID());
+
+ 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_u_op = _ip_data[ip].N_u_op;
+
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const& dNdx_p = _ip_data[ip].dNdx_p;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ auto& eps = _ip_data[ip].eps;
+ auto& S_L = _ip_data[ip].saturation;
+
+ 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 K_SR = _process_data.solid_bulk_modulus(t, x_position)[0];
+ auto const K_LR = _process_data.fluid_bulk_modulus(t, x_position)[0];
+ auto const temperature = _process_data.temperature(t, x_position)[0];
+ auto const porosity = _process_data.flow_material->getPorosity(
+ material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
+ auto const rho_LR = _process_data.flow_material->getFluidDensity(
+ -p_cap_ip, temperature);
+ auto const& b = _process_data.specific_body_force;
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value>::identity2;
+
+ S_L = _process_data.flow_material->getSaturation(
+ material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
+
+ double const dS_L_dp_cap =
+ _process_data.flow_material->getSaturationDerivative(
+ material_id, t, x_position, -p_cap_ip, temperature, S_L);
+
+ double const k_rel =
+ _process_data.flow_material->getRelativePermeability(
+ t, x_position, -p_cap_ip, temperature, S_L);
+ auto const mu = _process_data.flow_material->getFluidViscosity(
+ -p_cap_ip, temperature);
+
+ double const K_intrinsic =
+ _process_data.intrinsic_permeability(t, x_position)[0];
+ double const K_over_mu = K_intrinsic / mu;
+
+ //
+ // displacement equation, displacement part
+ //
+ eps.noalias() = B * u;
+
+ auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ temperature);
+
+ //
+ // displacement equation, displacement part
+ //
+ K.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_LR;
+ rhs.template segment<displacement_size>(displacement_index).noalias() +=
+ N_u_op.transpose() * rho * b * w;
+
+ //
+ // pressure equation, pressure part.
+ //
+ double const a0 = S_L * (alpha - porosity) / K_SR;
+ // Volumetric average specific storage of the solid and fluid phases.
+ double const specific_storage =
+ dS_L_dp_cap * (p_cap_ip * a0 - porosity) +
+ S_L * (porosity / K_LR + a0);
+ M.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
+
+ K.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() +=
+ dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+
+ rhs.template segment<pressure_size>(pressure_index).noalias() +=
+ dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+
+ //
+ // displacement equation, pressure part
+ //
+ K.template block<displacement_size, pressure_size>(displacement_index,
+ pressure_index)
+ .noalias() -= B.transpose() * alpha * S_L * identity2 * N_p * w;
+
+ //
+ // pressure equation, displacement part.
+ //
+ M.template block<pressure_size, displacement_size>(pressure_index,
+ displacement_index)
+ .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
+ identity2.transpose() * B * w;
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ 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)
+{
+ assert(local_x.size() == pressure_size + displacement_size);
+
+ auto p_L =
+ 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_L_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 =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypePressure::NodalMatrixType storage_p =
+ ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
+ pressure_size);
+
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size, pressure_size>
+ Kup = ShapeMatricesTypeDisplacement::template MatrixType<
+ displacement_size, pressure_size>::Zero(displacement_size,
+ pressure_size);
+
+ typename ShapeMatricesTypeDisplacement::template MatrixType<
+ pressure_size, displacement_size>
+ Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
+ pressure_size, displacement_size>::Zero(pressure_size,
+ displacement_size);
+
+ double const& dt = _process_data.dt;
+ auto const material_id =
+ _process_data.flow_material->getMaterialID(_element.getID());
+
+ 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_u_op = _ip_data[ip].N_u_op;
+
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const& dNdx_p = _ip_data[ip].dNdx_p;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
+
+ typename ShapeMatricesTypeDisplacement::GlobalDimVectorType u_ip(
+ DisplacementDim);
+ for (int i = 0; i < u_ip.size(); ++i)
+ {
+ NumLib::shapeFunctionInterpolate(
+ u.segment(i * ShapeFunctionDisplacement::NPOINTS,
+ ShapeFunctionDisplacement::NPOINTS),
+ N_u, u_ip.coeffRef(i));
+ }
+
+ auto& eps = _ip_data[ip].eps;
+ auto& S_L = _ip_data[ip].saturation;
+ auto const& sigma_eff = _ip_data[ip].sigma_eff;
+
+ 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 K_SR = _process_data.solid_bulk_modulus(t, x_position)[0];
+ auto const K_LR = _process_data.fluid_bulk_modulus(t, x_position)[0];
+ auto const temperature = _process_data.temperature(t, x_position)[0];
+
+ auto const porosity = _process_data.flow_material->getPorosity(
+ material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
+ auto const rho_LR = _process_data.flow_material->getFluidDensity(
+ -p_cap_ip, temperature);
+ auto const& b = _process_data.specific_body_force;
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value>::identity2;
+
+ S_L = _process_data.flow_material->getSaturation(
+ material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
+
+ double const dS_L_dp_cap =
+ _process_data.flow_material->getSaturationDerivative(
+ material_id, t, x_position, -p_cap_ip, temperature, S_L);
+
+ double const d2S_L_dp_cap_2 =
+ _process_data.flow_material->getSaturationDerivative2(
+ material_id, t, x_position, -p_cap_ip, temperature, S_L);
+
+ double const k_rel =
+ _process_data.flow_material->getRelativePermeability(
+ t, x_position, -p_cap_ip, temperature, S_L);
+ auto const mu = _process_data.flow_material->getFluidViscosity(
+ -p_cap_ip, temperature);
+
+ double const K_intrinsic =
+ _process_data.intrinsic_permeability(t, x_position)[0];
+ double const K_over_mu = K_intrinsic / mu;
+ //
+ // displacement equation, displacement part
+ //
+ eps.noalias() = B * u;
+
+ auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ temperature);
+
+ 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_LR;
+ local_rhs.template segment<displacement_size>(displacement_index)
+ .noalias() -=
+ (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+
+ //
+ // displacement equation, pressure part
+ //
+ Kup.noalias() += B.transpose() * alpha * S_L * identity2 * N_p * w;
+
+ /* For future implementation including swelling.
+ double const dsigma_eff_dp_cap = -K_intrinsic * m_swell * n *
+ std::pow(S_L, n - 1) * dS_L_dp_cap *
+ identity2;
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() -= B.transpose() * dsigma_eff_dp_cap * N_p * w;
+ */
+
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() -= B.transpose() * alpha *
+ (S_L + p_cap_ip * dS_L_dp_cap) * identity2 * N_p * w;
+
+ local_Jac
+ .template block<displacement_size, pressure_size>(
+ displacement_index, pressure_index)
+ .noalias() +=
+ N_u_op.transpose() * porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
+ //
+ // pressure equation, displacement part.
+ //
+ Kpu.noalias() += N_p.transpose() * S_L * rho_LR * alpha *
+ identity2.transpose() * B * w;
+
+ //
+ // pressure equation, pressure part.
+ //
+ laplace_p.noalias() +=
+ dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+
+ double const a0 = (alpha - porosity) / K_SR;
+ double const specific_storage =
+ dS_L_dp_cap * (p_cap_ip * S_L * a0 - porosity) +
+ S_L * (porosity / K_LR + S_L * a0);
+
+ double const dspecific_storage_dp_cap =
+ d2S_L_dp_cap_2 * (p_cap_ip * S_L * a0 - porosity) +
+ dS_L_dp_cap *
+ (porosity / K_LR + a0 * 3 * S_L + dS_L_dp_cap * p_cap_ip * a0);
+
+ storage_p.noalias() +=
+ N_p.transpose() * rho_LR * specific_storage * N_p * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N_p.transpose() * rho_LR * p_cap_dot_ip *
+ dspecific_storage_dp_cap * N_p * w;
+
+ /* In the derivation there is a div(du/dt) term in the Jacobian, but
+ * this implementation increases the total runtime by 1%. Maybe a very
+ * large step is needed to see the increase of efficiency.
+ double div_u_dot = 0;
+ for (int i = 0; i < DisplacementDim; ++i)
+ {
+ div_u_dot +=
+ (dNdx_u *
+ u_dot.template segment<ShapeFunctionDisplacement::NPOINTS>(
+ i * ShapeFunctionDisplacement::NPOINTS))[i];
+ }
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap * alpha *
+ div_u_dot * N_p * w;
+ */
+
+ double const dk_rel_dS_l =
+ _process_data.flow_material->getRelativePermeabilityDerivative(
+ t, x_position, -p_cap_ip, temperature, S_L);
+ typename ShapeMatricesTypeDisplacement::GlobalDimVectorType const
+ grad_p_cap = -dNdx_p * p_L;
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx_p.transpose() * rho_LR * K_over_mu * grad_p_cap *
+ dk_rel_dS_l * dS_L_dp_cap * N_p * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx_p.transpose() * rho_LR * rho_LR * K_over_mu * b *
+ dk_rel_dS_l * dS_L_dp_cap * N_p * w;
+
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
+ dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+ }
+
+ // 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() = Kpu / dt;
+
+ // pressure equation
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
+ laplace_p * p_L + storage_p * p_L_dot + Kpu * u_dot;
+
+ // displacement equation
+ local_rhs.template segment<displacement_size>(displacement_index)
+ .noalias() += Kup * p_L;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& RichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtSigma(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ static const int kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
+
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ auto const& sigma = _ip_data[ip].sigma_eff;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma);
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& RichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtEpsilon(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ auto const num_intpts = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, kelvin_vector_size, num_intpts);
+
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ auto const& eps = _ip_data[ip].eps;
+ cache_mat.col(ip) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+ }
+
+ return cache;
+}
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& RichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::getIntPtDarcyVelocity(const double t,
+ GlobalVector const&
+ current_solution,
+ NumLib::LocalToGlobalIndexMap const&
+ dof_table,
+ std::vector<double>& cache) const
+{
+ auto const num_intpts = _ip_data.size();
+
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_x = current_solution.get(indices);
+
+ cache.clear();
+ auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
+ double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, DisplacementDim, num_intpts);
+
+ auto p_L =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ pressure_size> const>(local_x.data() + pressure_index,
+ pressure_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_p = _ip_data[ip].N_p;
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+
+ auto const temperature = _process_data.temperature(t, x_position)[0];
+ auto const mu = _process_data.flow_material->getFluidViscosity(
+ -p_cap_ip, temperature);
+ double const K_over_mu =
+ _process_data.intrinsic_permeability(t, x_position)[0] / mu;
+ auto const rho_LR = _process_data.flow_material->getFluidDensity(
+ -p_cap_ip, temperature);
+ auto const& b = _process_data.specific_body_force;
+
+ // Compute the velocity
+ auto const& dNdx_p = _ip_data[ip].dNdx_p;
+ cache_matrix.col(ip).noalias() =
+ -K_over_mu * dNdx_p * p_L - K_over_mu * rho_LR * b;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+std::vector<double> const& RichardsMechanicsLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ getIntPtSaturation(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ auto const num_intpts = _ip_data.size();
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(cache, 1,
+ num_intpts);
+
+ for (unsigned ip = 0; ip < num_intpts; ++ip)
+ {
+ cache_mat[ip] = _ip_data[ip].saturation;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianForPressureEquations(
+ const double /*t*/, const std::vector<double>& /*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*/,
+ const LocalCoupledSolutions& /*local_coupled_solutions*/)
+{
+ OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianForDeformationEquations(
+ const double /*t*/, const std::vector<double>& /*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*/,
+ const LocalCoupledSolutions& /*local_coupled_solutions*/)
+{
+ OGS_FATAL("RichardsMechanics; The staggered scheme is not implemented.");
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobianForStaggeredScheme(
+ const double t,
+ const std::vector<double>& 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,
+ const LocalCoupledSolutions& local_coupled_solutions)
+{
+ // For the equations with pressure
+ if (local_coupled_solutions.process_id == 0)
+ {
+ assembleWithJacobianForPressureEquations(
+ t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+ local_b_data, local_Jac_data, local_coupled_solutions);
+ return;
+ }
+
+ // For the equations with deformation
+ assembleWithJacobianForDeformationEquations(
+ t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
+ local_b_data, local_Jac_data, local_coupled_solutions);
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ postNonLinearSolverConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const use_monolithic_scheme)
+{
+ const int displacement_offset =
+ use_monolithic_scheme ? displacement_index : 0;
+
+ auto u =
+ Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_offset,
+ displacement_size);
+ double const& dt = _process_data.dt;
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ int const n_integration_points = _integration_method.getNumberOfPoints();
+ for (int ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_u = _ip_data[ip].N_u;
+ auto const& dNdx_u = _ip_data[ip].dNdx_u;
+ auto const temperature = _process_data.temperature(t, x_position)[0];
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunctionDisplacement,
+ ShapeMatricesTypeDisplacement>(_element,
+ N_u);
+ auto const B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx_u, N_u, x_coord, _is_axially_symmetric);
+
+ auto& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ temperature);
+ }
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ typename IntegrationMethod, int DisplacementDim>
+void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
+ ShapeFunctionPressure, IntegrationMethod,
+ DisplacementDim>::
+ computeSecondaryVariableConcrete(double const t,
+ std::vector<double> const& local_x)
+{
+ auto p_L =
+ Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+ pressure_size> const>(local_x.data() + pressure_index,
+ pressure_size);
+ auto const material_id =
+ _process_data.flow_material->getMaterialID(_element.getID());
+
+ SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ double saturation_avg = 0;
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N_p = _ip_data[ip].N_p;
+ auto const temperature = _process_data.temperature(t, x_position)[0];
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+ auto& S_L = _ip_data[ip].saturation;
+ S_L = _process_data.flow_material->getSaturation(
+ material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
+ saturation_avg += S_L;
+ }
+ saturation_avg /= n_integration_points;
+
+ (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+}
+
+} // namespace RichardsMechanics
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
new file mode 100644
index 00000000000..ca91e763e88
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -0,0 +1,237 @@
+/**
+ * \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"
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
new file mode 100644
index 00000000000..d1ef1a3cd66
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -0,0 +1,381 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
new file mode 100644
index 00000000000..fe176bd3e3d
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.h
@@ -0,0 +1,140 @@
+/**
+ * \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
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
new file mode 100644
index 00000000000..df5a6edf8d0
--- /dev/null
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -0,0 +1,104 @@
+/**
+ * \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
--
GitLab