From fb19014ddffb552fce58b9794c46190861725dcc Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 27 Feb 2023 17:36:43 +0100
Subject: [PATCH] [PL/LD] LargeDeformation using MFront
Implementation of a finite strain mechanics process based solely on the
MFront material models.
---
Applications/ApplicationsLib/ProjectData.cpp | 32 ++
.../LARGE_DEFORMATION/c_LARGE_DEFORMATION.md | 1 +
.../LARGE_DEFORMATION/constitutive_relation | 1 +
.../process_variables/i_process_variables.md | 0
.../process_variables/t_process_variable.md | 1 +
.../LARGE_DEFORMATION/t_initial_stress.md | 1 +
.../t_reference_temperature.md | 1 +
.../t_specific_body_force.md | 1 +
ProcessLib/Deformation/GMatrixPolicy.h | 13 +-
ProcessLib/Deformation/NonLinearBMatrix.h | 120 +++++
ProcessLib/LargeDeformation/CMakeLists.txt | 16 +
.../ConstitutiveRelations/Base.h | 159 +++++++
.../ConstitutiveRelations/ConstitutiveData.h | 84 ++++
.../ConstitutiveModels.h | 38 ++
.../ConstitutiveSetting.cpp | 62 +++
.../ConstitutiveSetting.h | 42 ++
.../CreateConstitutiveSetting.cpp | 54 +++
.../CreateConstitutiveSetting.h | 31 ++
.../ConstitutiveRelations/Gravity.cpp | 23 +
.../ConstitutiveRelations/Gravity.h | 40 ++
.../ConstitutiveRelations/Invoke.h | 54 +++
.../ConstitutiveRelations/MaterialState.h | 32 ++
.../ConstitutiveRelations/SolidDensity.cpp | 26 ++
.../ConstitutiveRelations/SolidDensity.h | 24 +
.../ConstitutiveRelations/SolidMechanics.cpp | 85 ++++
.../ConstitutiveRelations/SolidMechanics.h | 85 ++++
.../ConstitutiveRelations/StressData.h | 32 ++
.../CreateLargeDeformationProcess.cpp | 162 +++++++
.../CreateLargeDeformationProcess.h | 85 ++++
.../LargeDeformation/CreateLocalAssemblers.h | 91 ++++
.../LargeDeformation/IntegrationPointData.h | 26 ++
.../LargeDeformation/LargeDeformationFEM.h | 415 ++++++++++++++++++
.../LargeDeformationProcess.cpp | 217 +++++++++
.../LargeDeformationProcess.h | 89 ++++
.../LargeDeformationProcessData.h | 62 +++
.../LocalAssemblerInterface.h | 176 ++++++++
scripts/cmake/ProcessesSetup.cmake | 4 +
37 files changed, 2379 insertions(+), 6 deletions(-)
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/c_LARGE_DEFORMATION.md
create mode 120000 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/constitutive_relation
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/i_process_variables.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/t_process_variable.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_initial_stress.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_reference_temperature.md
create mode 100644 Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_specific_body_force.md
create mode 100644 ProcessLib/Deformation/NonLinearBMatrix.h
create mode 100644 ProcessLib/LargeDeformation/CMakeLists.txt
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/Base.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveData.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveModels.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.cpp
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/Invoke.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/MaterialState.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.cpp
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.cpp
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.h
create mode 100644 ProcessLib/LargeDeformation/ConstitutiveRelations/StressData.h
create mode 100644 ProcessLib/LargeDeformation/CreateLargeDeformationProcess.cpp
create mode 100644 ProcessLib/LargeDeformation/CreateLargeDeformationProcess.h
create mode 100644 ProcessLib/LargeDeformation/CreateLocalAssemblers.h
create mode 100644 ProcessLib/LargeDeformation/IntegrationPointData.h
create mode 100644 ProcessLib/LargeDeformation/LargeDeformationFEM.h
create mode 100644 ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
create mode 100644 ProcessLib/LargeDeformation/LargeDeformationProcess.h
create mode 100644 ProcessLib/LargeDeformation/LargeDeformationProcessData.h
create mode 100644 ProcessLib/LargeDeformation/LocalAssemblerInterface.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index cf9658d5d8c..d1bbb257401 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -85,6 +85,9 @@
#ifdef OGS_BUILD_PROCESS_HYDROMECHANICS
#include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
#endif
+#ifdef OGS_BUILD_PROCESS_LARGEDEFORMATION
+#include "ProcessLib/LargeDeformation/CreateLargeDeformationProcess.h"
+#endif
#ifdef OGS_BUILD_PROCESS_LIE
#include "ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.h"
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
@@ -874,6 +877,35 @@ void ProjectData::parseProcesses(
}
else
#endif
+#ifdef OGS_BUILD_PROCESS_LARGEDEFORMATION
+ if (type == "LARGE_DEFORMATION")
+ {
+ switch (_mesh_vec[0]->getDimension())
+ {
+ case 2:
+ process = ProcessLib::LargeDeformation::
+ createLargeDeformationProcess<2>(
+ name, *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters,
+ _local_coordinate_system, integration_order,
+ process_config, _media);
+ break;
+ case 3:
+ process = ProcessLib::LargeDeformation::
+ createLargeDeformationProcess<3>(
+ name, *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters,
+ _local_coordinate_system, integration_order,
+ process_config, _media);
+ break;
+ default:
+ OGS_FATAL(
+ "LARGE_DEFORMATION process does not support given "
+ "dimension");
+ }
+ }
+ else
+#endif
#ifdef OGS_BUILD_PROCESS_LIE
if (type == "HYDRO_MECHANICS_WITH_LIE")
{
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/c_LARGE_DEFORMATION.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/c_LARGE_DEFORMATION.md
new file mode 100644
index 00000000000..5349d2c6274
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/c_LARGE_DEFORMATION.md
@@ -0,0 +1 @@
+Finite strain deformation process.
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/constitutive_relation b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/constitutive_relation
new file mode 120000
index 00000000000..cfcb5e03bbe
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/constitutive_relation
@@ -0,0 +1 @@
+../../../../material/solid/constitutive_relation
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/i_process_variables.md
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/t_process_variable.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/t_process_variable.md
new file mode 100644
index 00000000000..d55b23cfca7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/process_variables/t_process_variable.md
@@ -0,0 +1 @@
+The displacement vector field.
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_initial_stress.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_initial_stress.md
new file mode 100644
index 00000000000..9832658f886
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_initial_stress.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::SmallDeformation::SmallDeformationProcessData::initial_stress
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_reference_temperature.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_reference_temperature.md
new file mode 100644
index 00000000000..ce255b04b79
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_reference_temperature.md
@@ -0,0 +1 @@
+Reference temperature, an optional input.
diff --git a/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_specific_body_force.md
new file mode 100644
index 00000000000..a481e852d81
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/LARGE_DEFORMATION/t_specific_body_force.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::SmallDeformation::SmallDeformationProcessData::specific_body_force
diff --git a/ProcessLib/Deformation/GMatrixPolicy.h b/ProcessLib/Deformation/GMatrixPolicy.h
index 363222a461b..2665fd33521 100644
--- a/ProcessLib/Deformation/GMatrixPolicy.h
+++ b/ProcessLib/Deformation/GMatrixPolicy.h
@@ -20,11 +20,11 @@ template <typename ShapeFunction, int DisplacementDim>
class GMatrixPolicyType
{
private:
- /// Reusing the ShapeMatrixPolicy vector type.
+ /// Fixed size vector type independent of the ShapeMatrixPolicy needed for
+ /// storage of gradient vector in the MPL VariableArray.
template <int N>
- using VectorType =
- typename ShapeMatrixPolicyType<ShapeFunction,
- DisplacementDim>::template VectorType<N>;
+ using VectorTypeFixedSize = typename EigenFixedShapeMatrixPolicy<
+ ShapeFunction, DisplacementDim>::template VectorType<N>;
/// Reusing the ShapeMatrixPolicy matrix type.
template <int N, int M>
@@ -40,7 +40,8 @@ public:
using GradientMatrixType = MatrixType<DisplacementDim * DisplacementDim +
(DisplacementDim == 2 ? 1 : 0),
_number_of_dof>;
- using GradientVectorType = VectorType<DisplacementDim * DisplacementDim +
- (DisplacementDim == 2 ? 1 : 0)>;
+ using GradientVectorType =
+ VectorTypeFixedSize<DisplacementDim * DisplacementDim +
+ (DisplacementDim == 2 ? 1 : 0)>;
};
} // namespace ProcessLib
diff --git a/ProcessLib/Deformation/NonLinearBMatrix.h b/ProcessLib/Deformation/NonLinearBMatrix.h
new file mode 100644
index 00000000000..4b191a41849
--- /dev/null
+++ b/ProcessLib/Deformation/NonLinearBMatrix.h
@@ -0,0 +1,120 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <cmath>
+
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+
+namespace ProcessLib
+{
+namespace NonLinearBMatrix
+{
+/// Fills a non linear B-matrix based on given shape function dN/dx values and
+/// displacement gradient.
+template <int DisplacementDim,
+ int NPOINTS,
+ typename BMatrixType,
+ typename GradientVectorType,
+ typename N_Type,
+ typename DNDX_Type>
+BMatrixType computeBMatrix(DNDX_Type const& dNdx,
+ N_Type const& N,
+ GradientVectorType const& grad_u,
+ const double radius,
+ const bool is_axially_symmetric)
+{
+ static_assert(0 < DisplacementDim && DisplacementDim <= 3,
+ "NonLinearBMatrix::computeBMatrix: DisplacementDim must be "
+ "in range [1,3].");
+
+ BMatrixType B = BMatrixType::Zero(
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
+ NPOINTS * DisplacementDim);
+
+ switch (DisplacementDim)
+ {
+ case 3:
+ for (int i = 0; i < NPOINTS; ++i)
+ {
+ B(0, i) = dNdx(0, i) * grad_u[0];
+ B(1, i) = dNdx(1, i) * grad_u[1];
+ B(2, i) = dNdx(2, i) * grad_u[2];
+
+ B(0, NPOINTS + i) = dNdx(0, i) * grad_u[3];
+ B(1, NPOINTS + i) = dNdx(1, i) * grad_u[4];
+ B(2, NPOINTS + i) = dNdx(2, i) * grad_u[5];
+
+ B(0, 2 * NPOINTS + i) = dNdx(0, i) * grad_u[6];
+ B(1, 2 * NPOINTS + i) = dNdx(1, i) * grad_u[7];
+ B(2, 2 * NPOINTS + i) = dNdx(2, i) * grad_u[8];
+
+ B(3, i) = (dNdx(1, i) * grad_u[0] + dNdx(0, i) * grad_u[1]) /
+ std::sqrt(2);
+ B(4, i) = (dNdx(2, i) * grad_u[1] + dNdx(1, i) * grad_u[2]) /
+ std::sqrt(2);
+ B(5, i) = (dNdx(2, i) * grad_u[0] + dNdx(0, i) * grad_u[2]) /
+ std::sqrt(2);
+
+ B(3, NPOINTS + i) =
+ (dNdx(1, i) * grad_u[3] + dNdx(0, i) * grad_u[4]) /
+ std::sqrt(2);
+ B(4, NPOINTS + i) =
+ (dNdx(2, i) * grad_u[4] + dNdx(1, i) * grad_u[5]) /
+ std::sqrt(2);
+ B(5, NPOINTS + i) =
+ (dNdx(2, i) * grad_u[3] + dNdx(0, i) * grad_u[5]) /
+ std::sqrt(2);
+
+ B(3, 2 * NPOINTS + i) =
+ (dNdx(1, i) * grad_u[6] + dNdx(0, i) * grad_u[7]) /
+ std::sqrt(2);
+ B(4, 2 * NPOINTS + i) =
+ (dNdx(2, i) * grad_u[7] + dNdx(1, i) * grad_u[8]) /
+ std::sqrt(2);
+ B(5, 2 * NPOINTS + i) =
+ (dNdx(2, i) * grad_u[6] + dNdx(0, i) * grad_u[8]) /
+ std::sqrt(2);
+ }
+ break;
+ case 2:
+ for (int i = 0; i < NPOINTS; ++i)
+ {
+ B(0, i) = dNdx(0, i) * grad_u[0];
+ B(1, i) = dNdx(1, i) * grad_u[1];
+ // B(2, i) = 0;
+
+ B(0, NPOINTS + i) = dNdx(0, i) * grad_u[2];
+ B(1, NPOINTS + i) = dNdx(1, i) * grad_u[3];
+ // B(2, NPOINTS + i) = 0;
+
+ B(3, i) = (dNdx(1, i) * grad_u[0] + dNdx(0, i) * grad_u[1]) /
+ std::sqrt(2);
+ B(3, NPOINTS + i) =
+ (dNdx(1, i) * grad_u[2] + dNdx(0, i) * grad_u[3]) /
+ std::sqrt(2);
+ }
+ if (is_axially_symmetric)
+ {
+ for (int i = 0; i < NPOINTS; ++i)
+ {
+ B(2, i) = grad_u[4] * N[i] / radius;
+ }
+ }
+ break;
+ default:
+ break;
+ }
+
+ return B;
+}
+} // namespace NonLinearBMatrix
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/CMakeLists.txt b/ProcessLib/LargeDeformation/CMakeLists.txt
new file mode 100644
index 00000000000..f8a62ae04ab
--- /dev/null
+++ b/ProcessLib/LargeDeformation/CMakeLists.txt
@@ -0,0 +1,16 @@
+get_source_files(SOURCES)
+append_source_files(SOURCES ConstitutiveRelations)
+
+if (NOT OGS_USE_MFRONT)
+ message(FATAL_ERROR
+"LargeDeformations process without MFront library is not available. Use enable \
+MFront in the CMake settings via OGS_USE_MFRONT variable.")
+endif()
+
+ogs_add_library(LargeDeformation ${SOURCES})
+target_link_libraries(LargeDeformation PUBLIC ProcessLib PRIVATE ParameterLib)
+
+target_precompile_headers(LargeDeformation PRIVATE [["BaseLib/Error.h"]]
+ [["BaseLib/ConfigTree.h"]] [["BaseLib/Logging.h"]]
+ [["ProcessLib/Process.h"]] [["MaterialLib/MPL/Medium.h"]]
+ [["MaterialLib/MPL/Property.h"]] <Eigen/Core>)
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/Base.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/Base.h
new file mode 100644
index 00000000000..04493186ac8
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/Base.h
@@ -0,0 +1,159 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "BaseLib/StrongType.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/Tensor.h"
+#include "MathLib/KelvinVector.h"
+#include "ParameterLib/SpatialPosition.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::LargeDeformation
+{
+template <int DisplacementDim>
+using KelvinVector = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+
+template <int DisplacementDim>
+using KelvinMatrix = MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+
+template <int DisplacementDim>
+using GlobalDimVector = Eigen::Vector<double, DisplacementDim>;
+
+template <int DisplacementDim>
+using GlobalDimMatrix =
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim, Eigen::RowMajor>;
+
+/// Convenience alias for not a number.
+static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
+
+/// Used to set a Kelvin vector to all not-a-number.
+template <int DisplacementDim>
+constexpr KelvinVector<DisplacementDim> KVnan()
+{
+ return KelvinVector<DisplacementDim>::Constant(nan);
+}
+
+/// Used to set a Kelvin matrix to all not-a-number.
+template <int DisplacementDim>
+constexpr KelvinMatrix<DisplacementDim> KMnan()
+{
+ return KelvinMatrix<DisplacementDim>::Constant(nan);
+}
+
+/// Used to set a D dimensional vector to all not-a-number.
+template <int D>
+constexpr GlobalDimVector<D> DVnan()
+{
+ return GlobalDimVector<D>::Constant(nan);
+}
+
+/// Used to set a D x D matrix to all not-a-number.
+template <int D>
+constexpr GlobalDimMatrix<D> DMnan()
+{
+ return GlobalDimMatrix<D>::Constant(nan);
+}
+
+/// Used to set a Kelvin vector to all zero.
+template <int DisplacementDim>
+constexpr KelvinVector<DisplacementDim> KVzero()
+{
+ return KelvinVector<DisplacementDim>::Zero();
+}
+
+/// Used to set a Kelvin matrix to all zero.
+template <int DisplacementDim>
+constexpr KelvinMatrix<DisplacementDim> KMzero()
+{
+ return KelvinMatrix<DisplacementDim>::Zero();
+}
+
+/// Represents a previous state of type T.
+template <typename T>
+struct PrevState
+{
+ PrevState() = default;
+ explicit PrevState(T const& t) : t{t} {}
+ explicit PrevState(T&& t) : t{std::move(t)} {}
+
+ PrevState<T>& operator=(T const& u)
+ {
+ t = u;
+ return *this;
+ }
+
+ PrevState<T>& operator=(T&& u)
+ {
+ t = std::move(u);
+ return *this;
+ }
+
+ T& operator*() { return t; }
+ T const& operator*() const { return t; }
+
+ T* operator->() { return &t; }
+ T const* operator->() const { return &t; }
+
+private:
+ T t;
+};
+
+struct SpaceTimeData
+{
+ ParameterLib::SpatialPosition x;
+ double t;
+ double dt;
+};
+
+struct MediaData
+{
+ explicit MediaData(MaterialPropertyLib::Medium const& medium)
+ : medium{medium}, solid{medium.phase("Solid")}
+ {
+ }
+
+ MaterialPropertyLib::Medium const& medium;
+ MaterialPropertyLib::Phase const& solid;
+};
+
+using Temperature = BaseLib::StrongType<double, struct TemperatureTag>;
+
+template <int DisplacementDim>
+struct DeformationGradientData
+{
+ // TODO Move initialization to the local assembler.
+ MaterialPropertyLib::Tensor<DisplacementDim> deformation_gradient =
+ MaterialPropertyLib::Tensor<DisplacementDim>::Zero();
+
+ static auto reflect()
+ {
+ using Self = DeformationGradientData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName(
+ "deformation_gradient", &Self::deformation_gradient);
+ }
+};
+
+template <int DisplacementDim>
+struct StrainData
+{
+ // TODO Move initialization to the local assembler.
+ KelvinVector<DisplacementDim> eps = KVnan<DisplacementDim>();
+
+ static auto reflect()
+ {
+ using Self = StrainData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName("epsilon", &Self::eps);
+ }
+};
+
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveData.h
new file mode 100644
index 00000000000..e1a21bf1b4c
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveData.h
@@ -0,0 +1,84 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "Gravity.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
+#include "SolidDensity.h"
+#include "SolidMechanics.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+/// Data whose state must be tracked by the process.
+template <int DisplacementDim>
+struct StatefulData
+{
+ StressData<DisplacementDim> stress_data;
+
+ static auto reflect()
+ {
+ using Self = StatefulData<DisplacementDim>;
+
+ return Reflection::reflectWithoutName(&Self::stress_data);
+ }
+};
+
+/// Data whose state must be tracked by the process.
+template <int DisplacementDim>
+struct StatefulDataPrev
+{
+ PrevState<StressData<DisplacementDim>> stress_data;
+
+ StatefulDataPrev<DisplacementDim>& operator=(
+ StatefulData<DisplacementDim> const& state)
+ {
+ stress_data = state.stress_data;
+
+ return *this;
+ }
+};
+
+/// Data that is needed for output purposes, but not directly for the assembly.
+template <int DisplacementDim>
+struct OutputData
+{
+ StrainData<DisplacementDim> eps_data;
+
+ static auto reflect()
+ {
+ using Self = OutputData<DisplacementDim>;
+
+ return Reflection::reflectWithoutName(&Self::eps_data);
+ }
+};
+
+/// Data that is needed for the equation system assembly.
+template <int DisplacementDim>
+struct ConstitutiveData
+{
+ SolidMechanicsDataStateless<DisplacementDim> s_mech_data;
+ VolumetricBodyForce<DisplacementDim> volumetric_body_force;
+};
+
+/// Data that stores intermediate values, which are not needed outside the
+/// constitutive setting.
+template <int DisplacementDim>
+struct ConstitutiveTempData
+{
+ PrevState<StrainData<DisplacementDim>> eps_data_prev;
+ DeformationGradientData<DisplacementDim> deformation_gradient_data;
+ PrevState<DeformationGradientData<DisplacementDim>>
+ deformation_gradient_data_prev;
+ SolidDensity rho_SR;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveModels.h
new file mode 100644
index 00000000000..688a591c480
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveModels.h
@@ -0,0 +1,38 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "Gravity.h"
+#include "SolidDensity.h"
+#include "SolidMechanics.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+/// Constitutive models used for assembly.
+template <int DisplacementDim>
+struct ConstitutiveModels
+{
+ template <typename TRMProcessData>
+ explicit ConstitutiveModels(
+ TRMProcessData const& process_data,
+ SolidConstitutiveRelation<DisplacementDim> const& solid_material)
+ : s_mech_model(solid_material),
+ gravity_model(process_data.specific_body_force)
+ {
+ }
+
+ SolidMechanicsModel<DisplacementDim> s_mech_model;
+ SolidDensityModel rho_S_model;
+ GravityModel<DisplacementDim> gravity_model;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
new file mode 100644
index 00000000000..87b7cde859c
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.cpp
@@ -0,0 +1,62 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "ConstitutiveSetting.h"
+
+#include "Invoke.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void ConstitutiveSetting<DisplacementDim>::eval(
+ ConstitutiveModels<DisplacementDim>& models, double const t,
+ double const dt, ParameterLib::SpatialPosition const& x_position,
+ MaterialPropertyLib::Medium const& medium, double const T_ref,
+ GradientVectorType const& deformation_gradient,
+ GradientVectorType const& deformation_gradient_prev,
+ StatefulData<DisplacementDim>& state,
+ StatefulDataPrev<DisplacementDim> const& prev_state,
+ MaterialStateData<DisplacementDim>& mat_state,
+ ConstitutiveTempData<DisplacementDim>& tmp,
+ ConstitutiveData<DisplacementDim>& cd) const
+{
+ namespace MPL = MaterialPropertyLib;
+
+ auto& deformation_gradient_data = tmp.deformation_gradient_data;
+ deformation_gradient_data.deformation_gradient = deformation_gradient;
+ auto& deformation_gradient_data_prev = tmp.deformation_gradient_data_prev;
+ deformation_gradient_data_prev->deformation_gradient =
+ deformation_gradient_prev;
+ auto& rho_SR = tmp.rho_SR;
+
+ auto& s_mech_data = cd.s_mech_data;
+ auto& volumetric_body_force = cd.volumetric_body_force;
+
+ Temperature const T{T_ref};
+ SpaceTimeData const x_t{x_position, t, dt};
+ MediaData const media_data{medium};
+
+ assertEvalArgsUnique(models.s_mech_model);
+ models.s_mech_model.eval(
+ x_t, T, deformation_gradient_data, deformation_gradient_data_prev,
+ mat_state, prev_state.stress_data, state.stress_data, s_mech_data);
+
+ assertEvalArgsUnique(models.rho_S_model);
+ models.rho_S_model.eval(x_t, media_data, T, rho_SR);
+
+ assertEvalArgsUnique(models.gravity_model);
+ models.gravity_model.eval(rho_SR, volumetric_body_force);
+}
+
+template struct ConstitutiveSetting<2>;
+template struct ConstitutiveSetting<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.h
new file mode 100644
index 00000000000..63862f0167f
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/ConstitutiveSetting.h
@@ -0,0 +1,42 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "ConstitutiveData.h"
+#include "ConstitutiveModels.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct ConstitutiveSetting
+{
+ using GradientVectorType = Eigen::Matrix<
+ double,
+ DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0), 1>;
+
+ /// Evaluate the constitutive setting.
+ void eval(ConstitutiveModels<DisplacementDim>& models, double const t,
+ double const dt, ParameterLib::SpatialPosition const& x_position,
+ MaterialPropertyLib::Medium const& medium, double const T_ref,
+ GradientVectorType const& deformation_gradient,
+ GradientVectorType const& deformation_gradient_prev,
+ StatefulData<DisplacementDim>& state,
+ StatefulDataPrev<DisplacementDim> const& prev_state,
+ MaterialStateData<DisplacementDim>& mat_state,
+ ConstitutiveTempData<DisplacementDim>& tmp,
+ ConstitutiveData<DisplacementDim>& cd) const;
+};
+
+extern template struct ConstitutiveSetting<2>;
+extern template struct ConstitutiveSetting<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp b/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp
new file mode 100644
index 00000000000..0ec6c150565
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.cpp
@@ -0,0 +1,54 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "CreateConstitutiveSetting.h"
+
+#include "MaterialLib/SolidModels/CreateConstitutiveRelationsGeneric.h"
+#include "MaterialLib/SolidModels/MFront/CreateMFrontGeneric.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+std::unique_ptr<SolidConstitutiveRelation<DisplacementDim>> createMFrontGeneric(
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ BaseLib::ConfigTree const& config)
+{
+ namespace MSM = MaterialLib::Solids::MFront;
+ using namespace boost::mp11;
+
+ bool const library_path_is_relative_to_prj_file = true;
+
+ return MSM::createMFrontGeneric<
+ DisplacementDim, mp_list<MSM::DeformationGradient>,
+ mp_list<MSM::SecondPiolaKirchhoffStress>, mp_list<MSM::Temperature>>(
+ parameters, local_coordinate_system, config,
+ library_path_is_relative_to_prj_file);
+}
+
+template <int DisplacementDim>
+std::map<int, std::unique_ptr<SolidConstitutiveRelation<DisplacementDim>>>
+CreateConstitutiveSetting<DisplacementDim>::createSolidConstitutiveRelations(
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ BaseLib::ConfigTree const& config)
+{
+ return MaterialLib::Solids::createConstitutiveRelationsGeneric(
+ parameters, local_coordinate_system, config,
+ createMFrontGeneric<DisplacementDim>);
+}
+
+template struct CreateConstitutiveSetting<2>;
+template struct CreateConstitutiveSetting<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h
new file mode 100644
index 00000000000..300f1303d5a
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/CreateConstitutiveSetting.h
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "SolidMechanics.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct CreateConstitutiveSetting
+{
+ static std::map<int,
+ std::unique_ptr<SolidConstitutiveRelation<DisplacementDim>>>
+ createSolidConstitutiveRelations(
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ BaseLib::ConfigTree const& config);
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.cpp b/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.cpp
new file mode 100644
index 00000000000..e12252c257d
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.cpp
@@ -0,0 +1,23 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "Gravity.h"
+
+namespace ProcessLib::LargeDeformation
+{
+template <int DisplacementDim>
+void GravityModel<DisplacementDim>::eval(
+ SolidDensity const& rho_SR, VolumetricBodyForce<DisplacementDim>& out) const
+{
+ *out = *rho_SR * specific_body_force_;
+}
+
+template struct GravityModel<2>;
+template struct GravityModel<3>;
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.h
new file mode 100644
index 00000000000..223e31ff79b
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/Gravity.h
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "BaseLib/StrongType.h"
+#include "SolidDensity.h"
+
+namespace ProcessLib::LargeDeformation
+{
+template <int DisplacementDim>
+using VolumetricBodyForce =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>, struct GravityTag>;
+
+template <int DisplacementDim>
+struct GravityModel
+{
+ explicit GravityModel(
+ Eigen::Vector<double, DisplacementDim> const& specific_body_force)
+ : specific_body_force_(specific_body_force)
+ {
+ }
+
+ void eval(SolidDensity const& rho_SR,
+ VolumetricBodyForce<DisplacementDim>& out) const;
+
+private:
+ // TODO (naumov) Do we need to store this for each integration point?
+ Eigen::Vector<double, DisplacementDim> const specific_body_force_;
+};
+
+extern template struct GravityModel<2>;
+extern template struct GravityModel<3>;
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/Invoke.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/Invoke.h
new file mode 100644
index 00000000000..5fd51486065
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/Invoke.h
@@ -0,0 +1,54 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <boost/mp11.hpp>
+
+namespace ProcessLib::LargeDeformation
+{
+namespace detail
+{
+template <typename Result, typename Class, typename... Args>
+constexpr bool areEvalArgumentTypesUnique(Result (Class::*)(Args...))
+{
+ using namespace boost::mp11;
+ return mp_is_set<
+ mp_transform<std::remove_cvref_t, mp_list<Args...>>>::value;
+}
+
+template <typename Result, typename Class, typename... Args>
+constexpr bool areEvalArgumentTypesUnique(Result (Class::*)(Args...) const)
+{
+ using namespace boost::mp11;
+ return mp_is_set<
+ mp_transform<std::remove_cvref_t, mp_list<Args...>>>::value;
+}
+} // namespace detail
+
+/// Checks whether the argument types of the eval() method of the given type T
+/// are unique.
+///
+/// Argument types differing only in constness, reference or volatility are
+/// considered equal.
+template <typename T>
+constexpr bool areEvalArgumentTypesUnique()
+{
+ return detail::areEvalArgumentTypesUnique(&T::eval);
+}
+
+/// Statically asserts that the argument types of the passed Model's eval()
+/// method are unique.
+template <typename Model>
+constexpr void assertEvalArgsUnique(Model const&)
+{
+ static_assert(areEvalArgumentTypesUnique<std::remove_cvref_t<Model>>());
+}
+
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/MaterialState.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/MaterialState.h
new file mode 100644
index 00000000000..771a7e56313
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/MaterialState.h
@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "MaterialLib/SolidModels/MechanicsBase.h"
+
+namespace ProcessLib::LargeDeformation
+{
+template <int DisplacementDim>
+class MaterialStateData
+{
+ using MSV = typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables;
+
+public:
+ explicit MaterialStateData(std::unique_ptr<MSV>&& material_state_variables)
+ : material_state_variables(std::move(material_state_variables))
+ {
+ }
+
+ void pushBackState() { material_state_variables->pushBackState(); }
+
+ std::unique_ptr<MSV> material_state_variables;
+};
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.cpp
new file mode 100644
index 00000000000..71716ad4828
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.cpp
@@ -0,0 +1,26 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "SolidDensity.h"
+
+namespace ProcessLib::LargeDeformation
+{
+void SolidDensityModel::eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ Temperature const& temperature,
+ SolidDensity& out) const
+{
+ namespace MPL = MaterialPropertyLib;
+ MPL::VariableArray variables;
+ variables.temperature = *temperature;
+
+ *out = media_data.solid.property(MPL::PropertyType::density)
+ .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+}
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.h
new file mode 100644
index 00000000000..b98043d0986
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidDensity.h
@@ -0,0 +1,24 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "Base.h"
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::LargeDeformation
+{
+using SolidDensity = BaseLib::StrongType<double, struct SolidDensityTag>;
+
+struct SolidDensityModel
+{
+ void eval(SpaceTimeData const& x_t, MediaData const& media_data,
+ Temperature const& temperature, SolidDensity& out) const;
+};
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.cpp b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.cpp
new file mode 100644
index 00000000000..87b7cd5b90c
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.cpp
@@ -0,0 +1,85 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "SolidMechanics.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void SolidMechanicsModel<DisplacementDim>::eval(
+ SpaceTimeData const& x_t,
+ Temperature const& temperature,
+ DeformationGradientData<DisplacementDim> const& deformation_gradient_data,
+ PrevState<DeformationGradientData<DisplacementDim>> const&
+ deformation_gradient_data_prev,
+ MaterialStateData<DisplacementDim>& mat_state,
+ PrevState<StressData<DisplacementDim>> const& stress_data_prev,
+ StressData<DisplacementDim>& stress_data,
+ SolidMechanicsDataStateless<DisplacementDim>& current_stateless) const
+{
+ namespace MPL = MaterialPropertyLib;
+
+ // current state
+ MPL::VariableArray variables;
+ {
+ // thermodynamic forces
+ variables.stress = stress_data.sigma;
+
+ // gradient
+ variables.deformation_gradient =
+ deformation_gradient_data.deformation_gradient;
+
+ // external state variables
+ variables.temperature = *temperature;
+ }
+
+ // previous state
+ MPL::VariableArray variables_prev;
+ {
+ // thermodynamic forces
+ variables_prev.stress = stress_data_prev->sigma;
+
+ // gradient
+ variables_prev.deformation_gradient =
+ deformation_gradient_data_prev->deformation_gradient;
+
+ // external state variables
+ variables_prev.temperature = *temperature;
+ }
+
+ auto solution = solid_material_.integrateStress(
+ variables_prev, variables, x_t.t, x_t.x, x_t.dt,
+ *mat_state.material_state_variables);
+
+ if (!solution)
+ {
+ OGS_FATAL("Computation of local constitutive relation failed.");
+ }
+
+ auto& tdyn_forces_data = std::get<0>(*solution);
+
+ auto const view = solid_material_.createThermodynamicForcesView();
+
+ stress_data.sigma =
+ view.block(MSM::second_piola_kirchhoff_stress, tdyn_forces_data);
+ mat_state.material_state_variables = std::move(std::get<1>(*solution));
+
+ auto const& tangent_operator_data = std::get<2>(*solution);
+
+ current_stateless.stiffness_tensor = tangent_operator_blocks_view_.block(
+ MSM::second_piola_kirchhoff_stress, MSM::green_lagrange_strain,
+ tangent_operator_data);
+}
+
+template struct SolidMechanicsModel<2>;
+template struct SolidMechanicsModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.h
new file mode 100644
index 00000000000..1240a6687c7
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/SolidMechanics.h
@@ -0,0 +1,85 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "MaterialLib/SolidModels/MFront/MFrontGeneric.h"
+#include "MaterialLib/SolidModels/MFront/Variable.h"
+#include "MaterialState.h"
+#include "StressData.h"
+
+namespace ProcessLib::LargeDeformation
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct SolidMechanicsDataStateless
+{
+ KelvinMatrix<DisplacementDim> stiffness_tensor = KMnan<DisplacementDim>();
+ KelvinVector<DisplacementDim> J_uT_BT_K_N = KVnan<DisplacementDim>();
+ KelvinVector<DisplacementDim> J_up_BT_K_N = KVnan<DisplacementDim>();
+};
+
+namespace MSM = MaterialLib::Solids::MFront;
+
+template <int DisplacementDim>
+using SolidConstitutiveRelation =
+ MSM::MFrontGeneric<DisplacementDim,
+ boost::mp11::mp_list<MSM::DeformationGradient>,
+ boost::mp11::mp_list<MSM::SecondPiolaKirchhoffStress>,
+ boost::mp11::mp_list<MSM::Temperature>>;
+
+template <int DisplacementDim>
+struct SolidMechanicsModel
+{
+ explicit SolidMechanicsModel(
+ SolidConstitutiveRelation<DisplacementDim> const& solid_material)
+ : solid_material_(solid_material),
+ tangent_operator_blocks_view_{
+ solid_material.template createTangentOperatorBlocksView<
+ typename MSM::ForcesGradsCombinations<
+ boost::mp11::mp_list<MSM::GreenLagrangeStrain>,
+ boost::mp11::mp_list<MSM::SecondPiolaKirchhoffStress>,
+ boost::mp11::mp_list<MSM::Temperature>>::type>()}
+ {
+ }
+
+ void eval(
+ SpaceTimeData const& x_t,
+ Temperature const& temperature,
+ DeformationGradientData<DisplacementDim> const&
+ deformation_gradient_data,
+ PrevState<DeformationGradientData<DisplacementDim>> const&
+ deformation_gradient_data_prev,
+ MaterialStateData<DisplacementDim>& mat_state,
+ PrevState<StressData<DisplacementDim>> const& stress_data_prev,
+ StressData<DisplacementDim>& stress_data,
+ SolidMechanicsDataStateless<DisplacementDim>& current_stateless) const;
+
+ auto getInternalVariables() const
+ {
+ return solid_material_.getInternalVariables();
+ }
+
+private:
+ SolidConstitutiveRelation<DisplacementDim> const& solid_material_;
+
+ MSM::OGSMFrontTangentOperatorBlocksView<
+ DisplacementDim,
+ MSM::ForcesGradsCombinations<
+ boost::mp11::mp_list<MSM::GreenLagrangeStrain>,
+ boost::mp11::mp_list<MSM::SecondPiolaKirchhoffStress>,
+ boost::mp11::mp_list<MSM::Temperature>>::type>
+ tangent_operator_blocks_view_;
+};
+
+extern template struct SolidMechanicsModel<2>;
+extern template struct SolidMechanicsModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/ConstitutiveRelations/StressData.h b/ProcessLib/LargeDeformation/ConstitutiveRelations/StressData.h
new file mode 100644
index 00000000000..4007fd92350
--- /dev/null
+++ b/ProcessLib/LargeDeformation/ConstitutiveRelations/StressData.h
@@ -0,0 +1,32 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "Base.h"
+
+namespace ProcessLib::LargeDeformation
+{
+
+template <int DisplacementDim>
+struct StressData
+{
+ // Stress is stateful for some constitutive settings and therefore must be
+ // initialized to something valid, e.g., zero.
+ // TODO find a better solution for that.
+ KelvinVector<DisplacementDim> sigma = KVzero<DisplacementDim>();
+
+ static auto reflect()
+ {
+ using Self = StressData<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithName("sigma", &Self::sigma);
+ }
+};
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.cpp b/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.cpp
new file mode 100644
index 00000000000..6ef1479129c
--- /dev/null
+++ b/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.cpp
@@ -0,0 +1,162 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "CreateLargeDeformationProcess.h"
+
+#include <cassert>
+
+#include "ConstitutiveRelations/CreateConstitutiveSetting.h"
+#include "LargeDeformationProcess.h"
+#include "LargeDeformationProcessData.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
+#include "ParameterLib/Utils.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createLargeDeformationProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
+{
+ //! \ogs_file_param{prj__processes__process__type}
+ config.checkConfigParameter("type", "LARGE_DEFORMATION");
+ DBUG("Create LargeDeformationProcess.");
+
+ /// \section processvariablesld Process Variables
+
+ //! \ogs_file_param{prj__processes__process__LARGE_DEFORMATION__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ /// Primary process variables as they appear in the global component vector:
+ auto per_process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__LARGE_DEFORMATION__process_variables__process_variable}
+ "process_variable"});
+
+ DBUG("Associate displacement with process variable '{:s}'.",
+ per_process_variables.back().get().getName());
+
+ if (per_process_variables.back().get().getNumberOfGlobalComponents() !=
+ DisplacementDim)
+ {
+ OGS_FATAL(
+ "Number of components of the process variable '{:s}' is different "
+ "from the displacement dimension: got {:d}, expected {:d}",
+ per_process_variables.back().get().getName(),
+ per_process_variables.back().get().getNumberOfGlobalComponents(),
+ DisplacementDim);
+ }
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+ process_variables;
+ process_variables.push_back(std::move(per_process_variables));
+
+ /// \section parametersld Process Parameters
+ auto solid_constitutive_relations =
+ ConstitutiveRelations::CreateConstitutiveSetting<DisplacementDim>::
+ createSolidConstitutiveRelations(parameters,
+ local_coordinate_system, config);
+
+ // Specific body force
+ Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
+ {
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__LARGE_DEFORMATION__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 media_map =
+ MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+
+ // Reference temperature
+ auto const reference_temperature = ParameterLib::findOptionalTagParameter<
+ double>(
+ //! \ogs_file_param_special{prj__processes__process__LARGE_DEFORMATION__reference_temperature}
+ config, "reference_temperature", parameters, 1, &mesh);
+ if (reference_temperature)
+ {
+ DBUG("Use '{:s}' as reference temperature parameter.",
+ (*reference_temperature).name);
+ }
+
+ // Initial stress conditions
+ auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
+ //! \ogs_file_param_special{prj__processes__process__LARGE_DEFORMATION__initial_stress}
+ config, "initial_stress", parameters,
+ // Symmetric tensor size, 4 or 6, not a Kelvin vector.
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
+ &mesh);
+
+ LargeDeformationProcessData<DisplacementDim> process_data{
+ materialIDs(mesh),
+ std::move(media_map),
+ std::move(solid_constitutive_relations),
+ initial_stress,
+ specific_body_force,
+ reference_temperature};
+
+ SecondaryVariableCollection secondary_variables;
+
+ ProcessLib::createSecondaryVariables(config, secondary_variables);
+
+ return std::make_unique<LargeDeformationProcess<DisplacementDim>>(
+ std::move(name), mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(process_data), std::move(secondary_variables));
+}
+
+template std::unique_ptr<Process> createLargeDeformationProcess<2>(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+
+template std::unique_ptr<Process> createLargeDeformationProcess<3>(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.h b/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.h
new file mode 100644
index 00000000000..a911b6147a3
--- /dev/null
+++ b/ProcessLib/LargeDeformation/CreateLargeDeformationProcess.h
@@ -0,0 +1,85 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <optional>
+#include <string>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace ParameterLib
+{
+struct CoordinateSystem;
+struct ParameterBase;
+} // namespace ParameterLib
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+class Process;
+class ProcessVariable;
+} // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+std::unique_ptr<Process> createLargeDeformationProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+
+extern template std::unique_ptr<Process> createLargeDeformationProcess<2>(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+
+extern template std::unique_ptr<Process> createLargeDeformationProcess<3>(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ std::optional<ParameterLib::CoordinateSystem> const&
+ local_coordinate_system,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/CreateLocalAssemblers.h b/ProcessLib/LargeDeformation/CreateLocalAssemblers.h
new file mode 100644
index 00000000000..25df30ec75c
--- /dev/null
+++ b/ProcessLib/LargeDeformation/CreateLocalAssemblers.h
@@ -0,0 +1,91 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "BaseLib/Logging.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Utils/LocalAssemblerFactoryForDimGreaterEqualN.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+namespace detail
+{
+template <int GlobalDim,
+ template <typename /* shp fct */, int /* global dim */>
+ class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface,
+ IntegrationMethodProviderOrIntegrationOrder ProviderOrOrder,
+ typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ProviderOrOrder const& provider_or_order,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ DBUG("Create local assemblers.");
+
+ auto const& integration_method_provider =
+ getIntegrationMethodProvider(provider_or_order);
+
+ using IntMethProv =
+ std::remove_cvref_t<decltype(integration_method_provider)>;
+ using LocAsmFactory = ProcessLib::LocalAssemblerFactorySD<
+ LocalAssemblerInterface, LocalAssemblerImplementation, IntMethProv,
+ GlobalDim, ExtraCtorArgs...>;
+
+ LocAsmFactory factory(dof_table, integration_method_provider);
+ local_assemblers.resize(mesh_elements.size());
+
+ DBUG("Calling local assembler builder for all mesh elements.");
+ GlobalExecutor::transformDereferenced(
+ factory, 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 /* shp fct */, int /* global dim */>
+ class LocalAssemblerImplementation,
+ typename LocalAssemblerInterface,
+ IntegrationMethodProviderOrIntegrationOrder ProviderOrOrder,
+ typename... ExtraCtorArgs>
+void createLocalAssemblers(
+ std::vector<MeshLib::Element*> const& mesh_elements,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+ ProviderOrOrder const& provider_or_order,
+ ExtraCtorArgs&&... extra_ctor_args)
+{
+ DBUG("Create local assemblers.");
+
+ detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
+ dof_table, mesh_elements, local_assemblers, provider_or_order,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+} // namespace LargeDeformation
+
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/IntegrationPointData.h b/ProcessLib/LargeDeformation/IntegrationPointData.h
new file mode 100644
index 00000000000..ec0f0de8406
--- /dev/null
+++ b/ProcessLib/LargeDeformation/IntegrationPointData.h
@@ -0,0 +1,26 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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
+
+namespace ProcessLib::LargeDeformation
+{
+
+template <typename BMatricesType, typename ShapeMatricesType,
+ int DisplacementDim>
+struct IntegrationPointData final
+{
+ double integration_weight;
+ typename ShapeMatricesType::NodalRowVectorType N;
+ typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ProcessLib::LargeDeformation
diff --git a/ProcessLib/LargeDeformation/LargeDeformationFEM.h b/ProcessLib/LargeDeformation/LargeDeformationFEM.h
new file mode 100644
index 00000000000..34b68953855
--- /dev/null
+++ b/ProcessLib/LargeDeformation/LargeDeformationFEM.h
@@ -0,0 +1,415 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <limits>
+#include <memory>
+#include <vector>
+
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "MathLib/EigenBlockMatrixView.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/LocalDOF.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ParameterLib/Parameter.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/GMatrix.h"
+#include "ProcessLib/Deformation/GMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/Deformation/NonLinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+namespace MPL = MaterialPropertyLib;
+
+/// Used for the 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;
+};
+
+template <typename ShapeFunction, int DisplacementDim>
+class LargeDeformationLocalAssembler
+ : public LargeDeformationLocalAssemblerInterface<DisplacementDim>
+{
+public:
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+ using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;
+
+ using BMatrixType = typename BMatricesType::BMatrixType;
+ using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType;
+ using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
+ using NodalDisplacementVectorType =
+ typename BMatricesType::NodalForceVectorType;
+
+ using GMatricesType = GMatrixPolicyType<ShapeFunction, DisplacementDim>;
+ using GradientVectorType = typename GMatricesType::GradientVectorType;
+ using GradientMatrixType = typename GMatricesType::GradientMatrixType;
+ using IpData =
+ IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;
+
+ static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
+ DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>;
+ static constexpr auto& sigma_geom_op = MathLib::eigenBlockMatrixView<
+ DisplacementDim,
+ Eigen::Matrix<double, DisplacementDim, DisplacementDim>>;
+
+ LargeDeformationLocalAssembler(LargeDeformationLocalAssembler const&) =
+ delete;
+ LargeDeformationLocalAssembler(LargeDeformationLocalAssembler&&) = delete;
+
+ LargeDeformationLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ NumLib::GenericIntegrationMethod const& integration_method,
+ bool const is_axially_symmetric,
+ LargeDeformationProcessData<DisplacementDim>& process_data)
+ : LargeDeformationLocalAssemblerInterface<DisplacementDim>(
+ e, integration_method, is_axially_symmetric, process_data)
+ {
+ unsigned const n_integration_points =
+ this->integration_method_.getNumberOfPoints();
+
+ _ip_data.resize(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ DisplacementDim>(
+ e, is_axially_symmetric, this->integration_method_);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto& ip_data = _ip_data[ip];
+ auto const& sm = shape_matrices[ip];
+ _ip_data[ip].integration_weight =
+ this->integration_method_.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ;
+
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
+
+ _secondary_data.N[ip] = shape_matrices[ip].N;
+ }
+ }
+
+ void initializeConcrete() override
+ {
+ unsigned const n_integration_points =
+ this->integration_method_.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& ip_data = _ip_data[ip];
+
+ ParameterLib::SpatialPosition const x_position{
+ std::nullopt, this->element_.getID(), ip,
+ MathLib::Point3d(
+ NumLib::interpolateCoordinates<ShapeFunction,
+ ShapeMatricesType>(
+ this->element_, ip_data.N))};
+
+ /// Set initial stress from parameter.
+ if (this->process_data_.initial_stress != nullptr)
+ {
+ this->current_states_[ip].stress_data.sigma.noalias() =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector<
+ DisplacementDim>((*this->process_data_.initial_stress)(
+ std::numeric_limits<
+ double>::quiet_NaN() /* time independent */,
+ x_position));
+ }
+
+ double const t = 0; // TODO (naumov) pass t from top
+ auto& material_state = this->material_states_[ip];
+ this->solid_material_.initializeInternalStateVariables(
+ t, x_position, *material_state.material_state_variables);
+
+ material_state.pushBackState();
+ this->prev_states_[ip] = this->current_states_[ip];
+ }
+ }
+
+ typename ConstitutiveRelations::ConstitutiveData<DisplacementDim>
+ updateConstitutiveRelations(
+ Eigen::Ref<Eigen::VectorXd const> const& u,
+ Eigen::Ref<Eigen::VectorXd const> const& u_prev,
+ ParameterLib::SpatialPosition const& x_position, double const t,
+ double const dt,
+ IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>&
+ ip_data,
+ typename ConstitutiveRelations::ConstitutiveSetting<DisplacementDim>&
+ CS,
+ MaterialPropertyLib::Medium const& medium,
+ typename ConstitutiveRelations::StatefulData<DisplacementDim>&
+ current_state,
+ typename ConstitutiveRelations::StatefulDataPrev<DisplacementDim> const&
+ prev_state,
+ MaterialStateData<DisplacementDim>& material_state,
+ typename ConstitutiveRelations::OutputData<DisplacementDim>&
+ output_data) const
+ {
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+ auto const x_coord =
+ NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+ this->element_, N);
+
+ // For the 2D case the 33-component is needed (and the four entries
+ // of the non-symmetric matrix); In 3d there are nine entries.
+ GradientMatrixType G(
+ DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
+ ShapeFunction::NPOINTS * DisplacementDim);
+ Deformation::computeGMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
+ dNdx, G, this->is_axially_symmetric_, N, x_coord);
+
+ GradientVectorType const grad_u = G * u;
+ GradientVectorType gradient_vector_identity =
+ GradientVectorType::Zero();
+ if constexpr (DisplacementDim == 2)
+ {
+ gradient_vector_identity[0] = 1;
+ gradient_vector_identity[3] = 1;
+ gradient_vector_identity[4] = 1;
+ }
+ if constexpr (DisplacementDim == 3)
+ {
+ gradient_vector_identity[0] = 1;
+ gradient_vector_identity[4] = 1;
+ gradient_vector_identity[8] = 1;
+ }
+
+ auto const B_0 =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, x_coord, this->is_axially_symmetric_);
+
+ auto const B_N = NonLinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(dNdx, N, grad_u, x_coord,
+ this->is_axially_symmetric_);
+
+ auto const B = B_0 + 0.5 * B_N; // For Green-Lagrange strain.
+
+ double const T_ref =
+ this->process_data_.reference_temperature
+ ? (*this->process_data_.reference_temperature)(t, x_position)[0]
+ : std::numeric_limits<double>::quiet_NaN();
+
+ typename ConstitutiveRelations::ConstitutiveModels<DisplacementDim>
+ models(this->process_data_, this->solid_material_);
+ typename ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>
+ tmp;
+ typename ConstitutiveRelations::ConstitutiveData<DisplacementDim> CD;
+
+ output_data.eps_data.eps = B * u;
+
+ CS.eval(models, t, dt, x_position, //
+ medium, //
+ T_ref, //
+ grad_u + gradient_vector_identity, //
+ G * u_prev + gradient_vector_identity, //
+ current_state, prev_state, material_state, tmp, CD);
+
+ return CD;
+ }
+
+ void assemble(double const /*t*/, double const /*dt*/,
+ std::vector<double> const& /*local_x*/,
+ std::vector<double> const& /*local_x_prev*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& /*local_b_data*/) override
+ {
+ OGS_FATAL(
+ "LargeDeformationLocalAssembler: assembly without jacobian is not "
+ "implemented.");
+ }
+
+ void assembleWithJacobian(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_b_data,
+ std::vector<double>& local_Jac_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+
+ auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
+ local_b_data, local_matrix_size);
+
+ auto [u] = localDOF(local_x);
+ auto [u_prev] = localDOF(local_x_prev);
+
+ unsigned const n_integration_points =
+ this->integration_method_.getNumberOfPoints();
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(this->element_.getID());
+
+ typename ConstitutiveRelations::ConstitutiveSetting<DisplacementDim>
+ constitutive_setting;
+ auto const& medium =
+ *this->process_data_.media_map.getMedium(this->element_.getID());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ NumLib::interpolateXCoordinate<ShapeFunction,
+ ShapeMatricesType>(
+ this->element_, N);
+
+ // For the 2D case the 33-component is needed (and the four entries
+ // of the non-symmetric matrix); In 3d there are nine entries.
+ GradientMatrixType G(DisplacementDim * DisplacementDim +
+ (DisplacementDim == 2 ? 1 : 0),
+ ShapeFunction::NPOINTS * DisplacementDim);
+ Deformation::computeGMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS>(
+ dNdx, G, this->is_axially_symmetric_, N, x_coord);
+
+ GradientVectorType const grad_u = G * u;
+
+ auto const B_0 = LinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, x_coord, this->is_axially_symmetric_);
+
+ auto const B_N = NonLinearBMatrix::computeBMatrix<
+ DisplacementDim, ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, grad_u, x_coord, this->is_axially_symmetric_);
+
+ auto const B = B_0 + B_N;
+
+ auto const CD = updateConstitutiveRelations(
+ u, u_prev, x_position, t, dt, _ip_data[ip],
+ constitutive_setting, medium, this->current_states_[ip],
+ this->prev_states_[ip], this->material_states_[ip],
+ this->output_data_[ip]);
+
+ auto const& sigma = this->current_states_[ip].stress_data.sigma;
+ auto const& b = *CD.volumetric_body_force;
+ auto const& C = CD.s_mech_data.stiffness_tensor;
+
+ local_b.noalias() -=
+ (B.transpose() * sigma - N_u_op(N).transpose() * b) * w;
+
+ auto const sigma_tensor =
+ MathLib::KelvinVector::kelvinVectorToTensor(sigma);
+ if constexpr (DisplacementDim == 2)
+ {
+ using SigmaGeom =
+ typename ShapeMatricesType::template MatrixType<5, 5>;
+ SigmaGeom sigma_geom = SigmaGeom::Zero(5, 5);
+ sigma_geom.template block<4, 4>(0, 0) = sigma_geom_op(
+ sigma_tensor.template block<2, 2>(0, 0).eval());
+ sigma_geom(4, 4) = sigma_tensor(2, 2);
+
+ local_Jac.noalias() += G.transpose() * sigma_geom * G * w;
+ }
+ if constexpr (DisplacementDim == 3)
+ {
+ local_Jac.noalias() +=
+ G.transpose() * sigma_geom_op(sigma_tensor) * G * w;
+ }
+ local_Jac.noalias() += B.transpose() * C * B * w;
+ }
+ }
+
+ void postTimestepConcrete(Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_prev,
+ double const t, double const dt,
+ bool const /*use_monolithic_scheme*/,
+ int const /*process_id*/) override
+ {
+ unsigned const n_integration_points =
+ this->integration_method_.getNumberOfPoints();
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(this->element_.getID());
+
+ typename ConstitutiveRelations::ConstitutiveSetting<DisplacementDim>
+ constitutive_setting;
+
+ auto& medium =
+ *this->process_data_.media_map.getMedium(this->element_.getID());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ updateConstitutiveRelations(
+ local_x, local_x_prev, x_position, t, dt, _ip_data[ip],
+ constitutive_setting, medium, this->current_states_[ip],
+ this->prev_states_[ip], this->material_states_[ip],
+ this->output_data_[ip]);
+
+ this->material_states_[ip].pushBackState();
+ }
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ this->prev_states_[ip] = this->current_states_[ip];
+ }
+ }
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _secondary_data.N[integration_point];
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+private:
+ static constexpr auto localDOF(std::vector<double> const& x)
+ {
+ return NumLib::localDOF<
+ NumLib::Vectorial<ShapeFunction, DisplacementDim>>(x);
+ }
+
+private:
+ std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+
+ SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
+
+ static const int displacement_size =
+ ShapeFunction::NPOINTS * DisplacementDim;
+};
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
new file mode 100644
index 00000000000..02733399248
--- /dev/null
+++ b/ProcessLib/LargeDeformation/LargeDeformationProcess.cpp
@@ -0,0 +1,217 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "LargeDeformationProcess.h"
+
+#include <cassert>
+
+#include "CreateLocalAssemblers.h"
+#include "LargeDeformationFEM.h"
+#include "MeshLib/Utils/IntegrationPointWriter.h"
+#include "MeshLib/Utils/getOrCreateMeshProperty.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/Reflection/ReflectionForExtrapolation.h"
+#include "ProcessLib/Reflection/ReflectionForIPWriters.h"
+#include "ProcessLib/Utils/SetIPDataInitialConditions.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+LargeDeformationProcess<DisplacementDim>::LargeDeformationProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ LargeDeformationProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables)
+ : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables)),
+ _process_data(std::move(process_data))
+{
+ _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+ _process_data.principal_stress_vector[0] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_1",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_vector[1] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_2",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_vector[2] =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_3",
+ MeshLib::MeshItemType::Cell, 3);
+
+ _process_data.principal_stress_values =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "principal_stress_values",
+ MeshLib::MeshItemType::Cell, 3);
+
+ ProcessLib::Reflection::addReflectedIntegrationPointWriters<
+ DisplacementDim>(LargeDeformationLocalAssemblerInterface<
+ DisplacementDim>::getReflectionDataForOutput(),
+ _integration_point_writer, integration_order,
+ _local_assemblers);
+}
+
+template <int DisplacementDim>
+bool LargeDeformationProcess<DisplacementDim>::isLinear() const
+{
+ return false;
+}
+
+template <int DisplacementDim>
+void LargeDeformationProcess<DisplacementDim>::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::LargeDeformation::createLocalAssemblers<
+ DisplacementDim, LargeDeformationLocalAssembler>(
+ mesh.getElements(), dof_table, _local_assemblers,
+ NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
+ _process_data);
+
+ auto add_secondary_variable = [&](std::string const& name,
+ int const num_components,
+ auto get_ip_values_function)
+ {
+ _secondary_variables.addSecondaryVariable(
+ name,
+ makeExtrapolator(num_components, getExtrapolator(),
+ _local_assemblers,
+ std::move(get_ip_values_function)));
+ };
+
+ ProcessLib::Reflection::addReflectedSecondaryVariables<DisplacementDim>(
+ LargeDeformationLocalAssemblerInterface<
+ DisplacementDim>::getReflectionDataForOutput(),
+ _secondary_variables, getExtrapolator(), _local_assemblers);
+
+ //
+ // enable output of internal variables defined by material models
+ //
+ ProcessLib::Deformation::solidMaterialInternalToSecondaryVariables<
+ LocalAssemblerInterface>(_process_data.solid_materials,
+ add_secondary_variable);
+
+ ProcessLib::Deformation::
+ solidMaterialInternalVariablesToIntegrationPointWriter(
+ _process_data.solid_materials, _local_assemblers,
+ _integration_point_writer, integration_order);
+
+ bool const remove_name_suffix = true;
+ setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
+ _local_assemblers, remove_name_suffix);
+
+ // Initialize local assemblers after all variables have been set.
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::initialize, _local_assemblers,
+ *_local_to_global_index_map);
+}
+
+template <int DisplacementDim>
+void LargeDeformationProcess<DisplacementDim>::assembleConcreteProcess(
+ double const t, double const dt, std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+ DBUG("Assemble LargeDeformationProcess.");
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ 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::executeSelectedMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ pv.getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, M, K,
+ b);
+
+ _global_output(t, process_id, M, K, b);
+}
+
+template <int DisplacementDim>
+void LargeDeformationProcess<DisplacementDim>::
+ assembleWithJacobianConcreteProcess(
+ double const t, double const dt, std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac)
+{
+ DBUG("AssembleWithJacobian LargeDeformationProcess.");
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ 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::executeSelectedMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x,
+ x_prev, process_id, M, K, b, Jac);
+
+ transformVariableFromGlobalVector(b, 0, *_local_to_global_index_map,
+ *_nodal_forces, std::negate<double>());
+
+ _global_output(t, process_id, M, K, b, &Jac);
+}
+
+template <int DisplacementDim>
+void LargeDeformationProcess<DisplacementDim>::postTimestepConcreteProcess(
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, double const t, double const dt,
+ int const process_id)
+{
+ DBUG("PostTimestep LargeDeformationProcess.");
+ std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+ dof_tables.reserve(x.size());
+ std::generate_n(std::back_inserter(dof_tables), x.size(),
+ [&]() { return _local_to_global_index_map.get(); });
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &LocalAssemblerInterface::postTimestep, _local_assemblers,
+ pv.getActiveElementIDs(), dof_tables, x, x_prev, t, dt,
+ _use_monolithic_scheme, process_id);
+}
+
+template <int DisplacementDim>
+void LargeDeformationProcess<DisplacementDim>::computeSecondaryVariableConcrete(
+ double const t, double const dt, std::vector<GlobalVector*> const& x,
+ GlobalVector const& x_prev, int const process_id)
+{
+ DBUG("Compute the secondary variables for LargeDeformationProcess.");
+ std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+ dof_tables.reserve(x.size());
+ std::generate_n(std::back_inserter(dof_tables), x.size(),
+ [&]() { return _local_to_global_index_map.get(); });
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
+ pv.getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id);
+}
+template class LargeDeformationProcess<2>;
+template class LargeDeformationProcess<3>;
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcess.h b/ProcessLib/LargeDeformation/LargeDeformationProcess.h
new file mode 100644
index 00000000000..3568e2fcc0f
--- /dev/null
+++ b/ProcessLib/LargeDeformation/LargeDeformationProcess.h
@@ -0,0 +1,89 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "LargeDeformationProcessData.h"
+#include "LocalAssemblerInterface.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+class LargeDeformationProcess final : public Process
+{
+public:
+ LargeDeformationProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ LargeDeformationProcessData<DisplacementDim>&& process_data,
+ SecondaryVariableCollection&& secondary_variables);
+
+ //! \name ODESystem interface
+ //! @{
+ bool isLinear() const override;
+ //! @}
+
+private:
+ using LocalAssemblerInterface =
+ LargeDeformationLocalAssemblerInterface<DisplacementDim>;
+
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override;
+
+ void assembleConcreteProcess(double const t, double const dt,
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
+ int const process_id, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b) override;
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, double const dt, std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
+ GlobalMatrix& Jac) override;
+
+ void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
+ double const t, double const dt,
+ int const process_id) override;
+
+ void computeSecondaryVariableConcrete(double const t, double const dt,
+ std::vector<GlobalVector*> const& x,
+ GlobalVector const& x_prev,
+ int const process_id) override;
+
+private:
+ LargeDeformationProcessData<DisplacementDim> _process_data;
+
+ std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+ MeshLib::PropertyVector<double>* _nodal_forces = nullptr;
+ MeshLib::PropertyVector<double>* _material_forces = nullptr;
+
+ Assembly::GlobalMatrixOutput _global_output;
+};
+
+extern template class LargeDeformationProcess<2>;
+extern template class LargeDeformationProcess<3>;
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/LargeDeformationProcessData.h b/ProcessLib/LargeDeformation/LargeDeformationProcessData.h
new file mode 100644
index 00000000000..45075f87305
--- /dev/null
+++ b/ProcessLib/LargeDeformation/LargeDeformationProcessData.h
@@ -0,0 +1,62 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <Eigen/Core>
+#include <memory>
+#include <utility>
+
+#include "ConstitutiveRelations/SolidMechanics.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+}
+} // namespace MaterialLib
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+struct LargeDeformationProcessData
+{
+ MeshLib::PropertyVector<int> const* const material_ids = nullptr;
+
+ MaterialPropertyLib::MaterialSpatialDistributionMap media_map;
+
+ std::map<int,
+ std::unique_ptr<ConstitutiveRelations::SolidConstitutiveRelation<
+ DisplacementDim>>>
+ solid_materials;
+
+ /// Optional, initial stress field. A symmetric tensor, short vector
+ /// representation of length 4 or 6, ParameterLib::Parameter<double>.
+ ParameterLib::Parameter<double> const* const initial_stress;
+
+ /// Specific body forces applied to the solid.
+ /// It is usually used to apply gravitational forces.
+ /// A vector of displacement dimension's length.
+ Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
+
+ ParameterLib::Parameter<double> const* const reference_temperature;
+
+ std::array<MeshLib::PropertyVector<double>*, 3> principal_stress_vector = {
+ nullptr, nullptr, nullptr};
+ MeshLib::PropertyVector<double>* principal_stress_values = nullptr;
+};
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/ProcessLib/LargeDeformation/LocalAssemblerInterface.h b/ProcessLib/LargeDeformation/LocalAssemblerInterface.h
new file mode 100644
index 00000000000..b8943377678
--- /dev/null
+++ b/ProcessLib/LargeDeformation/LocalAssemblerInterface.h
@@ -0,0 +1,176 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, 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 "ConstitutiveRelations/ConstitutiveData.h"
+#include "ConstitutiveRelations/ConstitutiveSetting.h"
+#include "ConstitutiveRelations/MaterialState.h"
+#include "LargeDeformationProcessData.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/Reflection/ReflectionSetIPData.h"
+#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+
+namespace ProcessLib
+{
+namespace LargeDeformation
+{
+template <int DisplacementDim>
+struct LargeDeformationLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+ LargeDeformationLocalAssemblerInterface(
+ MeshLib::Element const& e,
+ NumLib::GenericIntegrationMethod const& integration_method,
+ bool const is_axially_symmetric,
+ LargeDeformationProcessData<DisplacementDim>& process_data)
+ : process_data_(process_data),
+ integration_method_(integration_method),
+ element_(e),
+ is_axially_symmetric_(is_axially_symmetric),
+ solid_material_(MaterialLib::Solids::selectSolidConstitutiveRelation(
+ process_data_.solid_materials, process_data_.material_ids,
+ element_.getID()))
+ {
+ unsigned const n_integration_points =
+ integration_method_.getNumberOfPoints();
+
+ material_states_.reserve(n_integration_points);
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ material_states_.emplace_back(
+ solid_material_.createMaterialStateVariables());
+ }
+
+ current_states_.resize(n_integration_points);
+ prev_states_.resize(n_integration_points);
+ output_data_.resize(n_integration_points);
+ }
+ /// Returns number of read integration points.
+ std::size_t setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order)
+ {
+ if (integration_order !=
+ static_cast<int>(integration_method_.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element {:d} is different "
+ "from the integration order in the initial condition.",
+ element_.getID());
+ }
+
+ if (name.starts_with("material_state_variable_"))
+ {
+ std::string const variable_name = name.substr(24, name.size() - 24);
+ DBUG("Setting material state variable '{:s}'", variable_name);
+
+ auto const& internal_variables =
+ solid_material_.getInternalVariables();
+ if (auto const iv = std::find_if(
+ begin(internal_variables), end(internal_variables),
+ [&variable_name](auto const& iv)
+ { return iv.name == variable_name; });
+ iv != end(internal_variables))
+ {
+ return ProcessLib::
+ setIntegrationPointDataMaterialStateVariables(
+ values, material_states_,
+ &MaterialStateData<
+ DisplacementDim>::material_state_variables,
+ iv->reference);
+ }
+
+ WARN(
+ "Could not find variable {:s} in solid material model's "
+ "internal variables.",
+ variable_name);
+ return 0;
+ }
+
+ // TODO this logic could be pulled out of the local assembler into the
+ // process. That might lead to a slightly better performance due to less
+ // string comparisons.
+ return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
+ name, values, current_states_);
+ }
+
+ // TODO move to NumLib::ExtrapolatableElement
+ unsigned getNumberOfIntegrationPoints() const
+ {
+ return integration_method_.getNumberOfPoints();
+ }
+
+ int getMaterialID() const
+ {
+ return process_data_.material_ids == nullptr
+ ? 0
+ : (*process_data_.material_ids)[element_.getID()];
+ }
+
+ std::vector<double> getMaterialStateVariableInternalState(
+ std::function<std::span<double>(
+ typename MaterialLib::Solids::MechanicsBase<DisplacementDim>::
+ MaterialStateVariables&)> const& get_values_span,
+ int const& n_components) const
+ {
+ return ProcessLib::getIntegrationPointDataMaterialStateVariables(
+ material_states_,
+ &MaterialStateData<DisplacementDim>::material_state_variables,
+ get_values_span, n_components);
+ }
+
+ typename MaterialLib::Solids::MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const&
+ getMaterialStateVariablesAt(unsigned integration_point) const
+ {
+ return *material_states_[integration_point].material_state_variables;
+ }
+
+ static auto getReflectionDataForOutput()
+ {
+ using Self = LargeDeformationLocalAssemblerInterface<DisplacementDim>;
+
+ return ProcessLib::Reflection::reflectWithoutName(
+ &Self::current_states_, &Self::output_data_);
+ }
+
+protected:
+ LargeDeformationProcessData<DisplacementDim>& process_data_;
+
+ // Material state is special, because it contains both the current and the
+ // old state.
+ std::vector<MaterialStateData<DisplacementDim>> material_states_;
+ std::vector<typename ConstitutiveRelations::StatefulData<DisplacementDim>>
+ current_states_; // TODO maybe do not store but rather re-evaluate for
+ // state update
+ std::vector<
+ typename ConstitutiveRelations::StatefulDataPrev<DisplacementDim>>
+ prev_states_;
+ std::vector<typename ConstitutiveRelations::OutputData<DisplacementDim>>
+ output_data_;
+
+ NumLib::GenericIntegrationMethod const& integration_method_;
+ MeshLib::Element const& element_;
+ bool const is_axially_symmetric_;
+ typename ConstitutiveRelations::SolidConstitutiveRelation<
+ DisplacementDim> const& solid_material_;
+};
+
+} // namespace LargeDeformation
+} // namespace ProcessLib
diff --git a/scripts/cmake/ProcessesSetup.cmake b/scripts/cmake/ProcessesSetup.cmake
index 119f8e47c78..b80a7a86f64 100644
--- a/scripts/cmake/ProcessesSetup.cmake
+++ b/scripts/cmake/ProcessesSetup.cmake
@@ -26,6 +26,10 @@ set(_processes_list
TwoPhaseFlowWithPP
TwoPhaseFlowWithPrho
)
+if(OGS_USE_MFRONT)
+ set(_processes_list ${_processes_list} LargeDeformation)
+endif()
+
if(OGS_USE_PETSC)
set(_processes_list ${_processes_list} PhaseField)
endif()
--
GitLab