diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index d837a53b200fb4c6871649c8e36944ff3f78c9a9..d35f6cee634c14bf9981dbe01d92dd83a4c1bd0b 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -486,6 +486,47 @@ if(NOT OGS_USE_MPI)
cube_1e3_expected_pcs_0_ts_101_t_1.000000.vtu cube_1e3_pcs_0_ts_101_t_1.000000.vtu displacement displacement
)
+ # Mechanics; Small deformations, Ehlers (SDE)
+ AddTest(
+ NAME Mechanics_SDE_cube_1e0
+ PATH Mechanics/Ehlers
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e0.prj
+ TESTER vtkdiff
+ ABSTOL 1e-15 RELTOL 1e-13
+ DIFF_DATA
+ expected_cube_1e0_pcs_0_ts_101_t_2.550000.vtu cube_1e0_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+ expected_cube_1e0_pcs_0_ts_101_t_2.550000.vtu cube_1e0_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+ expected_cube_1e0_pcs_0_ts_203_t_5.100000.vtu cube_1e0_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+ expected_cube_1e0_pcs_0_ts_203_t_5.100000.vtu cube_1e0_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+ )
+ AddTest(
+ NAME Mechanics_SDE_cube_1e1
+ PATH Mechanics/Ehlers
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e1.prj
+ TESTER vtkdiff
+ ABSTOL 1e-15 RELTOL 1e-13
+ DIFF_DATA
+ expected_cube_1e1_pcs_0_ts_101_t_2.550000.vtu cube_1e1_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+ expected_cube_1e1_pcs_0_ts_101_t_2.550000.vtu cube_1e1_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+ expected_cube_1e1_pcs_0_ts_203_t_5.100000.vtu cube_1e1_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+ expected_cube_1e1_pcs_0_ts_203_t_5.100000.vtu cube_1e1_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+ )
+ AddTest(
+ NAME LARGE_Mechanics_SDE_cube_1e3
+ PATH Mechanics/Ehlers
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e3.prj
+ TESTER vtkdiff
+ ABSTOL 1e-15 RELTOL 1e-13
+ DIFF_DATA
+ expected_cube_1e3_pcs_0_ts_101_t_2.550000.vtu cube_1e3_pcs_0_ts_101_t_2.550000.vtu displacement displacement
+ expected_cube_1e3_pcs_0_ts_101_t_2.550000.vtu cube_1e3_pcs_0_ts_101_t_2.550000.vtu sigma_xx sigma_xx
+ expected_cube_1e3_pcs_0_ts_203_t_5.100000.vtu cube_1e3_pcs_0_ts_203_t_5.100000.vtu displacement displacement
+ expected_cube_1e3_pcs_0_ts_203_t_5.100000.vtu cube_1e3_pcs_0_ts_203_t_5.100000.vtu sigma_xx sigma_xx
+ )
+
# SQUARE 1x1 GROUNDWATER FLOW TEST -- AXIALLY SYMMETRIC
# test results are compared to 3D simulation on a wedge-shaped domain
AddTest(
diff --git a/Documentation/bibliography.bib b/Documentation/bibliography.bib
index a25118eb8f773d5d1950ea2e1f2a3937e4ec907f..a6e2af3e046e5208b219be9a7347b46df5535808 100644
--- a/Documentation/bibliography.bib
+++ b/Documentation/bibliography.bib
@@ -2,6 +2,17 @@
% Encoding: UTF-8
+@Article{Ehlers1995,
+ Title={A single-surface yield function for geomaterials},
+ Author={Ehlers, W},
+ Journal={Archive of Applied Mechanics},
+ Volume={65},
+ Number={4},
+ Pages={246--259},
+ Year={1995},
+ Publisher={Springer}
+}
+
@Book{Ericson:2004:RCD:1121584,
Title = {Real-Time Collision Detection (The Morgan Kaufmann Series in Interactive 3-D Technology) (The Morgan Kaufmann Series in Interactive 3D Technology)},
Author = {Ericson, Christer},
diff --git a/MaterialLib/SolidModels/CreateEhlers.h b/MaterialLib/SolidModels/CreateEhlers.h
new file mode 100644
index 0000000000000000000000000000000000000000..a03a55b0135d9f912a8cc81ae569eab2e013aef8
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateEhlers.h
@@ -0,0 +1,138 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_CREATEEHLERS_H_
+#define MATERIALLIB_SOLIDMODELS_CREATEEHLERS_H_
+
+#include <logog/include/logog.hpp>
+
+#include "MechanicsBase.h"
+#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+#include "Ehlers.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+std::unique_ptr<MechanicsBase<DisplacementDim>> createEhlers(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfigParameter("type", "Ehlers");
+ DBUG("Create Ehlers material");
+
+ auto& shear_modulus = ProcessLib::findParameter<double>(
+ config, "shear_modulus", parameters, 1);
+
+ DBUG("Use \'%s\' as shear modulus parameter.", shear_modulus.name.c_str());
+
+ auto& bulk_modulus = ProcessLib::findParameter<double>(
+ config, "bulk_modulus", parameters, 1);
+
+ DBUG("Use \'%s\' as bulk modulus parameter.", bulk_modulus.name.c_str());
+
+ auto& kappa = ProcessLib::findParameter<double>(
+ config, "kappa", parameters, 1);
+
+ DBUG("Use \'%s\' as kappa.",
+ kappa.name.c_str());
+
+ auto& beta = ProcessLib::findParameter<double>(
+ config, "beta", parameters, 1);
+
+ DBUG("Use \'%s\' as beta.",
+ beta.name.c_str());
+
+ auto& gamma = ProcessLib::findParameter<double>(
+ config, "gamma", parameters, 1);
+
+ DBUG("Use \'%s\' as gamma.",
+ gamma.name.c_str());
+
+ auto& hardening_modulus = ProcessLib::findParameter<double>(
+ config, "hardening_modulus", parameters, 1);
+
+ DBUG("Use \'%s\' as hardening modulus parameter.",
+ hardening_modulus.name.c_str());
+
+ auto& alpha = ProcessLib::findParameter<double>(
+ config, "alpha", parameters, 1);
+
+ DBUG("Use \'%s\' as alpha.",
+ alpha.name.c_str());
+
+ auto& delta = ProcessLib::findParameter<double>(
+ config, "delta", parameters, 1);
+
+ DBUG("Use \'%s\' as delta.",
+ delta.name.c_str());
+
+ auto& eps = ProcessLib::findParameter<double>(
+ config, "eps", parameters, 1);
+
+ DBUG("Use \'%s\' as eps.",
+ eps.name.c_str());
+
+ auto& m = ProcessLib::findParameter<double>(
+ config, "m", parameters, 1);
+
+ DBUG("Use \'%s\' as m.",
+ m.name.c_str());
+
+ auto& alphap = ProcessLib::findParameter<double>(
+ config, "alphap", parameters, 1);
+
+ DBUG("Use \'%s\' as alphap.",
+ alphap.name.c_str());
+
+ auto& deltap = ProcessLib::findParameter<double>(
+ config, "deltap", parameters, 1);
+
+ DBUG("Use \'%s\' as deltap.",
+ deltap.name.c_str());
+
+ auto& epsp = ProcessLib::findParameter<double>(
+ config, "epsp", parameters, 1);
+
+ DBUG("Use \'%s\' as epsp.",
+ epsp.name.c_str());
+
+ auto& paremeter_mp = ProcessLib::findParameter<double>(
+ config, "mp", parameters, 1);
+
+ DBUG("Use \'%s\' as mp.",
+ paremeter_mp.name.c_str());
+
+ auto& betap = ProcessLib::findParameter<double>(
+ config, "betap", parameters, 1);
+
+ DBUG("Use \'%s\' as betap.",
+ betap.name.c_str());
+
+ auto& gammap = ProcessLib::findParameter<double>(
+ config, "gammap", parameters, 1);
+
+ DBUG("Use \'%s\' as gammap.",
+ gammap.name.c_str());
+
+ typename SolidEhlers<DisplacementDim>::MaterialProperties mp{
+ shear_modulus, bulk_modulus, alpha, beta,
+ gamma, delta, eps, m,
+ alphap, betap, gammap, deltap,
+ epsp, paremeter_mp, kappa, hardening_modulus};
+
+ return std::unique_ptr<MechanicsBase<DisplacementDim>>{
+ new SolidEhlers<DisplacementDim>{mp}};
+}
+
+} // namespace Solids
+} // namespace MaterialLib
+
+#endif // MATERIALLIB_SOLIDMODELS_CREATEEHLERS_H_
diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..15a6e5fbeebc5ba99e1a120e81af0b77825c14f0
--- /dev/null
+++ b/MaterialLib/SolidModels/Ehlers-impl.h
@@ -0,0 +1,666 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+/**
+ * Common convenitions for naming:
+ * x_D - deviatoric part of tensor x
+ * x_V - volumetric part of tensor x
+ * x_p - a variable related to plastic potential
+ * x_prev - value of x in previous time step
+ *
+ * Variables used in the code:
+ * eps_D - deviatoric strain
+ * eps_p_D_dot - deviatoric increment of plastic strain
+ * eps_p_eff_dot - increment of effective plastic strain
+ * eps_p_V_dot - volumetric increment of plastic strain
+ * sigma_D_inverse_D - deviatoric part of sigma_D_inverse
+ *
+ * derivation of the flow rule
+ * theta - J3 / J2^(3 / 2) from yield function
+ * dtheta_dsigma - derivative of theta
+ * sqrtPhi - square root of Phi from plastic potential
+ * flow_D - deviatoric part of flow
+ * flow_V - volumetric part of flow
+ * lambda_flow_D - deviatoric increment of plastic strain
+ *
+ */
+#ifndef MATERIALLIB_SOLIDMODELS_EHLERS_IMPL_H_
+#define MATERIALLIB_SOLIDMODELS_EHLERS_IMPL_H_
+
+#include <boost/math/special_functions/pow.hpp>
+#include <logog/include/logog.hpp>
+#include "MaterialLib/SolidModels/KelvinVector.h"
+#include "NewtonRaphson.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct PhysicalStressWithInvariants final
+{
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+
+ explicit PhysicalStressWithInvariants(KelvinVector const& stress)
+ : value{stress},
+ D{Invariants::deviatoric_projection * stress},
+ I_1{Invariants::trace(stress)},
+ J_2{Invariants::J2(D)},
+ J_3{Invariants::J3(D)}
+ {
+ }
+
+ PhysicalStressWithInvariants(PhysicalStressWithInvariants const&) = default;
+ PhysicalStressWithInvariants& operator=(
+ PhysicalStressWithInvariants const&) = default;
+#if defined(_MSC_VER) && (_MSC_VER >= 1900)
+ PhysicalStressWithInvariants(PhysicalStressWithInvariants&&) = default;
+ PhysicalStressWithInvariants& operator=(PhysicalStressWithInvariants&&) =
+ default;
+#endif // _MSC_VER
+
+ KelvinVector value;
+ KelvinVector D;
+ double I_1;
+ double J_2;
+ double J_3;
+};
+
+/// Holds powers of 1 + gamma_p*theta to base 0, m_p, and m_p-1.
+struct OnePlusGamma_pTheta final
+{
+ OnePlusGamma_pTheta(double const gamma_p, double const theta,
+ double const m_p)
+ : value{1 + gamma_p * theta},
+ pow_m_p{std::pow(value, m_p)},
+ pow_m_p1{pow_m_p/value}
+ {
+ }
+
+ double const value;
+ double const pow_m_p;
+ double const pow_m_p1;
+};
+
+template <int DisplacementDim>
+double plasticFlowVolumetricPart(
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ double const sqrtPhi, double const alpha_p, double const beta_p,
+ double const delta_p, double const epsilon_p)
+{
+ return 3 * (alpha_p * s.I_1 +
+ 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) /
+ (2 * sqrtPhi) +
+ 3 * beta_p + 6 * epsilon_p * s.I_1;
+}
+
+template <int DisplacementDim>
+typename SolidEhlers<DisplacementDim>::KelvinVector
+plasticFlowDeviatoricPart(
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ OnePlusGamma_pTheta const& one_gt, double const sqrtPhi,
+ typename SolidEhlers<DisplacementDim>::KelvinVector const&
+ dtheta_dsigma,
+ double const gamma_p, double const m_p)
+{
+ return (one_gt.pow_m_p * (s.D +
+ s.J_2 * m_p * gamma_p * dtheta_dsigma / one_gt.value)) /
+ (2 * sqrtPhi);
+}
+template <int DisplacementDim>
+double yieldFunction(
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ typename SolidEhlers<DisplacementDim>::MaterialProperties const& mp,
+ double const t, ProcessLib::SpatialPosition const& x)
+{
+ double const alpha = mp.alpha(t, x)[0];
+ double const beta = mp.beta(t, x)[0];
+ double const delta = mp.delta(t, x)[0];
+ double const epsilon = mp.epsilon(t, x)[0];
+ double const gamma = mp.gamma(t, x)[0];
+ double const m = mp.m(t, x)[0];
+
+ double const I_1_squared = boost::math::pow<2>(s.I_1);
+ assert(s.J_2 != 0);
+
+ return std::sqrt(
+ s.J_2 * std::pow(1 +
+ gamma * s.J_3 /
+ boost::math::pow<3>(std::sqrt(s.J_2)),
+ m) +
+ alpha / 2. * I_1_squared +
+ boost::math::pow<2>(delta) * boost::math::pow<2>(I_1_squared)) +
+ beta * s.I_1 + epsilon * I_1_squared - mp.k;
+}
+
+template <int DisplacementDim>
+void calculatePlasticResidual(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ ProcessLib::KelvinVectorType<DisplacementDim> const& eps_D,
+ double const eps_V,
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ ProcessLib::KelvinVectorType<DisplacementDim> const& eps_p_D,
+ ProcessLib::KelvinVectorType<DisplacementDim> const& eps_p_D_dot,
+ double const eps_p_V,
+ double const eps_p_V_dot,
+ double const eps_p_eff_dot,
+ double const lambda,
+ typename SolidEhlers<DisplacementDim>::MaterialProperties const& _mp,
+ typename SolidEhlers<DisplacementDim>::ResidualVectorType& residual)
+{
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& identity2 = Invariants::identity2;
+
+ double const G = _mp.G(t, x)[0];
+ double const K = _mp.K(t, x)[0];
+
+ double const theta = s.J_3 / boost::math::pow<3>(std::sqrt(s.J_2));
+
+ // calculate stress residual
+ residual.template segment<KelvinVectorSize>(0).noalias() =
+ s.value / G - 2 * (eps_D - eps_p_D) -
+ K / G * (eps_V - eps_p_V) * identity2;
+
+ // deviatoric plastic strain
+ double const alpha_p = _mp.alpha_p(t, x)[0];
+ double const delta_p = _mp.delta_p(t, x)[0];
+ double const gamma_p = _mp.gamma_p(t, x)[0];
+ double const m_p = _mp.m_p(t, x)[0];
+ KelvinVector const sigma_D_inverse_D =
+ P_dev * MaterialLib::SolidModels::inverse(s.D);
+ KelvinVector const dtheta_dsigma =
+ theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
+
+ OnePlusGamma_pTheta const one_gt{gamma_p, theta, m_p};
+ double const sqrtPhi = std::sqrt(
+ s.J_2 * one_gt.pow_m_p + alpha_p / 2. * boost::math::pow<2>(s.I_1) +
+ boost::math::pow<2>(delta_p) * boost::math::pow<4>(s.I_1));
+ KelvinVector const flow_D = plasticFlowDeviatoricPart(
+ s, one_gt, sqrtPhi, dtheta_dsigma, gamma_p, m_p);
+ KelvinVector const lambda_flow_D = lambda * flow_D;
+
+ residual.template segment<KelvinVectorSize>(KelvinVectorSize)
+ .noalias() = eps_p_D_dot - lambda_flow_D;
+
+ // plastic volume strain
+ {
+ double const beta_p = _mp.beta_p(t, x)[0];
+ double const epsilon_p = _mp.epsilon_p(t, x)[0];
+
+ double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
+ s, sqrtPhi, alpha_p, beta_p, delta_p, epsilon_p);
+ residual(2 * KelvinVectorSize, 0) = eps_p_V_dot - lambda * flow_V;
+ }
+
+ // evolution of plastic equivalent strain
+ residual(2 * KelvinVectorSize + 1) =
+ eps_p_eff_dot -
+ std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
+
+ // yield function (for plastic multiplier)
+ residual(2 * KelvinVectorSize + 2) =
+ yieldFunction<DisplacementDim>(s, _mp, t, x) / G;
+}
+
+/// Special product of v with itself: v \odot v
+/// v is given in Kelvin mapping.
+template <int DisplacementDim>
+ProcessLib::KelvinMatrixType<DisplacementDim> s_odot_s(
+ ProcessLib::KelvinVectorType<DisplacementDim> const& v)
+{
+ ProcessLib::KelvinMatrixType<DisplacementDim> result;
+
+ result(0, 0) = v(0) * v(0);
+ result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
+ result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
+ result(0, 3) = result(3, 0) = v(0) * v(3);
+ result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
+ result(0, 5) = result(5, 0) = v(0) * v(5);
+
+ result(1, 1) = v(1) * v(1);
+ result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
+ result(1, 3) = result(3, 1) = v(3) * v(1);
+ result(1, 4) = result(4, 1) = v(1) * v(4);
+ result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
+
+ result(2, 2) = v(2) * v(2);
+ result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
+ result(2, 4) = result(4, 2) = v(4) * v(2);
+ result(2, 5) = result(5, 2) = v(5) * v(2);
+
+ result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
+ result(3, 4) = result(4, 3) =
+ v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
+ result(3, 5) = result(5, 3) =
+ v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
+
+ result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
+ result(4, 5) = result(5, 4) =
+ v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
+
+ result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
+ return result;
+}
+
+template <int DisplacementDim>
+void calculatePlasticJacobian(
+ double const dt,
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ typename SolidEhlers<DisplacementDim>::JacobianMatrix& jacobian,
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ double const lambda,
+ typename SolidEhlers<DisplacementDim>::MaterialProperties const& _mp)
+{
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& identity2 = Invariants::identity2;
+
+ double const G = _mp.G(t, x)[0];
+ double const K = _mp.K(t, x)[0];
+ double const m = _mp.m(t, x)[0];
+ double const alpha = _mp.alpha(t, x)[0];
+ double const beta = _mp.beta(t, x)[0];
+ double const gamma = _mp.gamma(t, x)[0];
+ double const delta = _mp.delta(t, x)[0];
+
+ double const alpha_p = _mp.alpha_p(t, x)[0];
+ double const beta_p = _mp.beta_p(t, x)[0];
+ double const gamma_p = _mp.gamma_p(t, x)[0];
+ double const delta_p = _mp.delta_p(t, x)[0];
+ double const epsilon_p = _mp.epsilon_p(t, x)[0];
+ double const m_p = _mp.m_p(t, x)[0];
+
+ double const theta = s.J_3 / boost::math::pow<3>(std::sqrt(s.J_2));
+ OnePlusGamma_pTheta const one_gt{gamma_p, theta, m_p};
+
+ // inverse of deviatoric stress tensor
+ if (Invariants::determinant(s.D) == 0)
+ {
+ OGS_FATAL("Determinant is zero. Matrix is non-invertable.");
+ }
+ // inverse of sigma_D
+ KelvinVector const sigma_D_inverse = MaterialLib::SolidModels::inverse(s.D);
+ KelvinVector const sigma_D_inverse_D = P_dev * sigma_D_inverse;
+
+ KelvinVector const dtheta_dsigma =
+ theta * sigma_D_inverse_D - 3. / 2. * theta / s.J_2 * s.D;
+
+ // deviatoric flow
+ double const sqrtPhi = std::sqrt(
+ s.J_2 * one_gt.pow_m_p + alpha_p / 2. * boost::math::pow<2>(s.I_1) +
+ boost::math::pow<2>(delta_p) * boost::math::pow<4>(s.I_1));
+ KelvinVector const flow_D = plasticFlowDeviatoricPart(
+ s, one_gt, sqrtPhi, dtheta_dsigma, gamma_p, m_p);
+ KelvinVector const lambda_flow_D = lambda * flow_D;
+
+ jacobian.setZero();
+
+ // G_11
+ jacobian.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
+ .noalias() = KelvinMatrix::Identity();
+
+ // G_12
+ jacobian
+ .template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
+ .noalias() = 2 * KelvinMatrix::Identity();
+
+ // G_13
+ jacobian.template block<KelvinVectorSize, 1>(0, 2 * KelvinVectorSize)
+ .noalias() = K / G * identity2;
+
+ // G_14 and G_15 are zero
+
+ // G_21 -- derivative of deviatoric flow
+
+ double const gm_p = gamma_p * m_p;
+ // intermediate variable for derivative of deviatoric flow
+ KelvinVector const M0 = s.J_2 / one_gt.value * dtheta_dsigma;
+ // derivative of Phi w.r.t. sigma
+ KelvinVector const dPhi_dsigma =
+ one_gt.pow_m_p * (s.D + gm_p * M0) +
+ (alpha_p * s.I_1 +
+ 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) *
+ identity2;
+
+ // intermediate variable for derivative of deviatoric flow
+ KelvinMatrix const M1 =
+ one_gt.pow_m_p *
+ (s.D * dPhi_dsigma.transpose() + gm_p * M0 * dPhi_dsigma.transpose());
+ // intermediate variable for derivative of deviatoric flow
+ KelvinMatrix const M2 =
+ one_gt.pow_m_p * (P_dev + s.D * gm_p * M0.transpose());
+ // second derivative of theta
+ KelvinMatrix const d2theta_dsigma2 =
+ theta * P_dev * s_odot_s<DisplacementDim>(sigma_D_inverse) * P_dev +
+ sigma_D_inverse_D * dtheta_dsigma.transpose() -
+ 3. / 2. * theta / s.J_2 * P_dev -
+ 3. / 2. * dtheta_dsigma / s.J_2 * s.D.transpose() +
+ 3. / 2. * theta / boost::math::pow<2>(s.J_2) * s.D * s.D.transpose();
+
+ // intermediate variable for derivative of deviatoric flow
+ KelvinMatrix const M3 =
+ gm_p * one_gt.pow_m_p1 *
+ ((s.D + (gm_p - gamma_p) * M0) * dtheta_dsigma.transpose() +
+ s.J_2 * d2theta_dsigma2);
+
+ // derivative of flow_D w.r.t. sigma
+ KelvinMatrix const dflow_D_dsigma =
+ (-M1 / (4 * boost::math::pow<3>(sqrtPhi)) + (M2 + M3) / (2 * sqrtPhi)) *
+ G;
+ jacobian
+ .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
+ .noalias() = -lambda * dflow_D_dsigma;
+
+ // G_22
+ jacobian
+ .template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
+ KelvinVectorSize)
+ .noalias() = KelvinMatrix::Identity() / dt;
+
+ // G_23 and G_24 are zero
+
+ // G_25
+ jacobian
+ .template block<KelvinVectorSize, 1>(KelvinVectorSize,
+ 2 * KelvinVectorSize + 2)
+ .noalias() = -flow_D;
+
+ // G_31
+ {
+ // derivative of flow_V w.r.t. sigma
+ KelvinVector const dflow_V_dsigma =
+ 3 * G *
+ (-(alpha_p * s.I_1 +
+ 4 * boost::math::pow<2>(delta_p) * boost::math::pow<3>(s.I_1)) /
+ (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
+ (alpha_p * identity2 +
+ 12 * boost::math::pow<2>(delta_p * s.I_1) * identity2) /
+ (2 * sqrtPhi) +
+ 2 * epsilon_p * identity2);
+
+ jacobian.template block<1, KelvinVectorSize>(2 * KelvinVectorSize, 0)
+ .noalias() = -lambda * dflow_V_dsigma.transpose();
+ }
+
+ // G_32 is zero
+
+ // G_33
+ jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize) = 1. / dt;
+
+ // G_34 is zero
+
+ // G_35
+ {
+ double const flow_V = plasticFlowVolumetricPart<DisplacementDim>(
+ s, sqrtPhi, alpha_p, beta_p, delta_p, epsilon_p);
+ jacobian(2 * KelvinVectorSize, 2 * KelvinVectorSize + 2) = -flow_V;
+ }
+
+ // increment of effectiv plastic strain
+ double const eff_flow =
+ std::sqrt(2. / 3. * lambda_flow_D.transpose() * lambda_flow_D);
+
+ if (eff_flow > 0)
+ {
+ // intermediate variable for derivative of plastic jacobian
+ KelvinVector const eff_flow23_lambda_flow_D =
+ -2 / 3. / eff_flow * lambda_flow_D;
+ // G_41
+ jacobian
+ .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 1, 0)
+ .noalias() = lambda * dflow_D_dsigma * eff_flow23_lambda_flow_D;
+ // G_45
+ jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 2) =
+ eff_flow23_lambda_flow_D.transpose() * flow_D;
+ }
+
+ // G_42 and G_43 are zero
+
+ // G_44
+ jacobian(2 * KelvinVectorSize + 1, 2 * KelvinVectorSize + 1) = 1. / dt;
+
+ // G_51
+ {
+ double const one_gt_pow_m = std::pow(one_gt.value, m);
+ double const gm = gamma * m;
+ // derivative of yield function w.r.t. sigma
+ KelvinVector const dF_dsigma =
+ G * (one_gt_pow_m * (s.D + gm * M0) +
+ (alpha * s.I_1 +
+ 4 * boost::math::pow<2>(delta) * boost::math::pow<3>(s.I_1)) *
+ identity2) /
+ (2. * sqrtPhi) +
+ G * (beta + 2 * epsilon_p * s.I_1) * identity2;
+
+ jacobian
+ .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
+ .noalias() = dF_dsigma.transpose() / G;
+ }
+
+ // G_54
+ jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
+ -_mp.kappa(t, x)[0] * _mp.hardening_coefficient(t, x)[0] / G;
+
+ // G_52, G_53, G_55 are zero
+}
+
+/// Calculates the derivative of the residuals with respect to total
+/// strain. Implementation fully implicit only.
+template <int DisplacementDim>
+ProcessLib::KelvinMatrixType<DisplacementDim> calculateDResidualDEps(
+ double const K, double const G)
+{
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& P_sph = Invariants::spherical_projection;
+ auto const& I = ProcessLib::KelvinMatrixType<DisplacementDim>::Identity();
+
+ return -2. * I * P_dev - 3. * K / G * I * P_sph;
+}
+
+template <int DisplacementDim>
+void SolidEhlers<DisplacementDim>::MaterialProperties::
+ calculateIsotropicHardening(double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const eps_p_eff)
+{
+ k = kappa(t, x)[0] * (1. + eps_p_eff * hardening_coefficient(t, x)[0]);
+}
+
+template <int DisplacementDim>
+typename SolidEhlers<DisplacementDim>::KelvinVector predict_sigma(
+ double const G, double const K,
+ typename SolidEhlers<DisplacementDim>::KelvinVector const&
+ sigma_prev,
+ typename SolidEhlers<DisplacementDim>::KelvinVector const& eps,
+ typename SolidEhlers<DisplacementDim>::KelvinVector const& eps_prev,
+ double const eps_V)
+{
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+ auto const& P_dev = Invariants::deviatoric_projection;
+
+ // dimensionless initial hydrostatic pressure
+ double const pressure_prev = Invariants::trace(sigma_prev) / (-3. * G);
+ // initial strain invariant
+ double const e_prev = Invariants::trace(eps_prev);
+ // dimensioness hydrostatic stress increment
+ double const pressure = pressure_prev - K / G * (eps_V - e_prev);
+ // dimensionless deviatoric initial stress
+ typename SolidEhlers<DisplacementDim>::KelvinVector const
+ sigma_D_prev = P_dev * sigma_prev / G;
+ // dimensionless deviatoric stress
+ typename SolidEhlers<DisplacementDim>::KelvinVector const sigma_D =
+ sigma_D_prev + 2 * P_dev * (eps - eps_prev);
+ return sigma_D - pressure * Invariants::identity2;
+}
+
+template <int DisplacementDim>
+bool SolidEhlers<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ KelvinVector& sigma_final,
+ KelvinMatrix& C,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables)
+{
+ assert(dynamic_cast<MaterialStateVariables*>(&material_state_variables) !=
+ nullptr);
+ MaterialStateVariables& _state =
+ static_cast<MaterialStateVariables&>(material_state_variables);
+ _state.loadState();
+
+ using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
+ C.setZero();
+
+ // volumetric strain
+ double const eps_V = Invariants::trace(eps);
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ // deviatoric strain
+ KelvinVector const eps_D = P_dev * eps;
+
+ // dimensionless stress/hydrostatic pressure
+ double const G = _mp.G(t, x)[0];
+ double const K = _mp.K(t, x)[0];
+
+ KelvinVector sigma =
+ predict_sigma<DisplacementDim>(G, K, sigma_prev, eps, eps_prev, eps_V);
+
+ // update parameter
+ _mp.calculateIsotropicHardening(t, x, _state.eps_p_eff);
+
+ PhysicalStressWithInvariants<DisplacementDim> s{G * sigma};
+ // Quit early if sigma is zero (nothing to do) or if we are still in elastic
+ // zone.
+ if (sigma.squaredNorm() == 0 ||
+ yieldFunction<DisplacementDim>(s, _mp, t, x) < 0)
+ {
+ C.setZero();
+ C.template topLeftCorner<3, 3>().setConstant(K - 2. / 3 * G);
+ C.noalias() += 2 * G * KelvinMatrix::Identity();
+ }
+ else
+ {
+ JacobianMatrix jacobian;
+
+ // Linear solver for the newton loop is required after the loop with the
+ // same matrix. This saves one decomposition.
+ Eigen::FullPivLU<JacobianMatrix> linear_solver;
+
+ { // Newton loop for return mapping calculation.
+ auto const update_residual = [&](ResidualVectorType& residual) {
+
+ KelvinVector const eps_p_D_dot =
+ (_state.eps_p_D - _state.eps_p_D_prev) / dt;
+ double const eps_p_V_dot =
+ (_state.eps_p_V - _state.eps_p_V_prev) / dt;
+ double const eps_p_eff_dot =
+ (_state.eps_p_eff - _state.eps_p_eff_prev) / dt;
+ calculatePlasticResidual<DisplacementDim>(
+ t, x, eps_D, eps_V, s, _state.eps_p_D, eps_p_D_dot,
+ _state.eps_p_V, eps_p_V_dot, eps_p_eff_dot, _state.lambda,
+ _mp, residual);
+ };
+
+ auto const update_jacobian = [&](JacobianMatrix& jacobian) {
+ calculatePlasticJacobian<DisplacementDim>(dt, t, x, jacobian, s,
+ _state.lambda, _mp);
+ };
+
+ auto const update_solution = [&](
+ ResidualVectorType const& increment) {
+ sigma.noalias() +=
+ increment.template segment<KelvinVectorSize>(
+ KelvinVectorSize * 0);
+ s = PhysicalStressWithInvariants<DisplacementDim>{G * sigma};
+ _state.eps_p_D.noalias() +=
+ increment.template segment<KelvinVectorSize>(
+ KelvinVectorSize * 1);
+ _state.eps_p_V += increment(KelvinVectorSize * 2);
+ _state.eps_p_eff += increment(KelvinVectorSize * 2 + 1);
+ _state.lambda += increment(KelvinVectorSize * 2 + 2);
+
+ _mp.calculateIsotropicHardening(t, x, _state.eps_p_eff);
+ };
+
+ // TODO Make the following choice of maximum iterations and
+ // convergence criteria available from the input file configuration:
+ int const maximum_iterations(100);
+ double const tolerance(1e-14);
+
+ auto newton_solver =
+ NewtonRaphson<decltype(linear_solver), JacobianMatrix,
+ decltype(update_jacobian), ResidualVectorType,
+ decltype(update_residual),
+ decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual,
+ update_solution, maximum_iterations, tolerance);
+
+ auto const success_iterations = newton_solver.solve(jacobian);
+
+ if (!success_iterations)
+ return false;
+
+ // If the Newton loop didn't run, the linear solver will not be
+ // initialized.
+ // This happens usually for the first iteration of the first
+ // timestep.
+ if (*success_iterations == 0)
+ linear_solver.compute(jacobian);
+ }
+
+ // Calculate residual derivative w.r.t. strain
+ Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
+ Eigen::RowMajor>
+ dresidual_deps =
+ Eigen::Matrix<double, JacobianResidualSize, KelvinVectorSize,
+ Eigen::RowMajor>::Zero();
+ dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
+ .noalias() = calculateDResidualDEps<DisplacementDim>(K, G);
+
+ // Extract consistent tangent.
+ C.noalias() =
+ _mp.G(t, x)[0] *
+ linear_solver.solve(-dresidual_deps)
+ .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
+ }
+
+ // Update sigma.
+ sigma_final.noalias() = G * sigma;
+ return true;
+}
+
+} // namespace Solids
+} // namespace MaterialLib
+
+#endif // MATERIALLIB_SOLIDMODELS_EHLERS_IMPL_H_
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d494ce535768ad13b336451b5fea8368eaae5f92
--- /dev/null
+++ b/MaterialLib/SolidModels/Ehlers.cpp
@@ -0,0 +1,20 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "Ehlers.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template class SolidEhlers<2>;
+template class SolidEhlers<3>;
+
+} // namespace Solids
+} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/Ehlers.h b/MaterialLib/SolidModels/Ehlers.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6054564208f0f632f8145b6714fa4c392578b31
--- /dev/null
+++ b/MaterialLib/SolidModels/Ehlers.h
@@ -0,0 +1,217 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+/*
+ * Implementation of Ehler's single-surface model.
+ * see Ehler's paper "A single-surface yield function for geomaterials" for more
+ * details. \cite{Ehlers1995}
+ *
+ * Refer to "Single-surface benchmark of OpenGeoSys documentation
+ * (https://docs.opengeosys.org/docs/benchmarks/small-deformations/mechanics-plasticity-single-surface)"
+ * for more details for the tests.
+ */
+
+#ifndef MATERIALLIB_SOLIDMODELS_EHLERS_H_
+#define MATERIALLIB_SOLIDMODELS_EHLERS_H_
+
+#include <cfloat>
+#include <memory>
+#ifndef NDEBUG
+#include <ostream>
+#endif
+
+#include <Eigen/Dense>
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/Error.h"
+#include "KelvinVector.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+
+template <int DisplacementDim>
+class SolidEhlers final : public MechanicsBase<DisplacementDim>
+{
+public:
+ static int const KelvinVectorSize =
+ ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
+ static int const JacobianResidualSize =
+ 2 * KelvinVectorSize + 3; // 2 is the number of components in the
+ // jacobian/residual, not the space
+ // dimension. And 3 is for additional
+ // variables.
+ using ResidualVectorType = Eigen::Matrix<double, JacobianResidualSize, 1>;
+ using JacobianMatrix = Eigen::Matrix<double, JacobianResidualSize,
+ JacobianResidualSize, Eigen::RowMajor>;
+
+public:
+ //
+ // Variables specific to the material model.
+ //
+ struct MaterialProperties
+ {
+ using P = ProcessLib::Parameter<double>;
+
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+ /// material parameters in relation to Ehler's single-surface model
+ /// see Ehler's paper "A single-surface yield function for geomaterials"
+ /// for more details.
+ MaterialProperties(P const& G_, P const& K_, P const& alpha_,
+ P const& beta_, P const& gamma_, P const& delta_,
+ P const& epsilon_, P const& m_, P const& alpha_p_,
+ P const& beta_p_, P const& gamma_p_,
+ P const& delta_p_, P const& epsilon_p_,
+ P const& m_p_, P const& kappa_,
+ P const& hardening_coefficient_)
+ : G(G_),
+ K(K_),
+ alpha(alpha_),
+ beta(beta_),
+ gamma(gamma_),
+ delta(delta_),
+ epsilon(epsilon_),
+ m(m_),
+ alpha_p(alpha_p_),
+ beta_p(beta_p_),
+ gamma_p(gamma_p_),
+ delta_p(delta_p_),
+ epsilon_p(epsilon_p_),
+ m_p(m_p_),
+ kappa(kappa_),
+ hardening_coefficient(hardening_coefficient_)
+ {
+ }
+ // basic material parameters
+ P const& G; ///< shear modulus
+ P const& K; ///< bulk modulus
+
+ P const& alpha;
+ P const& beta;
+ P const& gamma;
+ P const& delta;
+ P const& epsilon;
+ P const& m;
+
+ P const& alpha_p;
+ P const& beta_p;
+ P const& gamma_p;
+ P const& delta_p;
+ P const& epsilon_p;
+ P const& m_p;
+
+ P const& kappa;
+ P const& hardening_coefficient;
+ // Drucker-Prager: Import kappa and beta in terms of Drucker-Prager
+ // criterion solution dependent values
+ double k;
+
+ void calculateIsotropicHardening(double const t,
+ ProcessLib::SpatialPosition const& x,
+ const double e_pv_curr);
+ };
+
+ struct MaterialStateVariables
+ : public MechanicsBase<DisplacementDim>::MaterialStateVariables
+ {
+ MaterialStateVariables()
+ : eps_p_D(KelvinVector::Zero()), eps_p_D_prev(KelvinVector::Zero())
+ {
+ }
+
+ void loadState()
+ {
+ eps_p_D = eps_p_D_prev;
+ eps_p_V = eps_p_V_prev;
+ eps_p_eff = eps_p_eff_prev;
+ lambda = 0;
+ }
+
+ void pushBackState() override
+ {
+ eps_p_D_prev = eps_p_D;
+ eps_p_V_prev = eps_p_V;
+ eps_p_eff_prev = eps_p_eff; // effective part of trace(eps_p)
+ lambda = 0;
+ }
+
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+
+ KelvinVector eps_p_D; ///< deviatoric plastic strain
+ double eps_p_V = 0; ///< volumetric strain
+ double eps_p_eff = 0; ///< effective plastic strain
+
+ // Initial values from previous timestep
+ KelvinVector eps_p_D_prev; ///< \copydoc eps_p_D
+ double eps_p_V_prev = 0; ///< \copydoc eps_p_V
+ double eps_p_eff_prev = 0; ///< \copydoc eps_p_eff
+
+ double lambda = 0; ///< plastic multiplier
+
+#ifndef NDEBUG
+ friend std::ostream& operator<<(std::ostream& os,
+ MaterialStateVariables const& m)
+ {
+ os << "State:\n"
+ << "eps_p_D: " << m.eps_p_D << "\n"
+ << "e_p: " << m.e_p << "\n"
+ << "eps_p_eff: " << m.eps_p_eff << "\n"
+ << "eps_p_D_prev: " << m.eps_p_D_prev << "\n"
+ << "e_p_prev: " << m.e_p_prev << "\n"
+ << "eps_p_eff_prev: " << m.eps_p_eff_prev << "\n"
+ << "lambda: " << m.lambda << "\n";
+ return os;
+ }
+#endif // NDEBUG
+ };
+
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() override
+ {
+ return std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+ new MaterialStateVariables};
+ }
+
+public:
+ using KelvinVector = ProcessLib::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix = ProcessLib::KelvinMatrixType<DisplacementDim>;
+
+public:
+ explicit SolidEhlers(MaterialProperties const& material_properties)
+ : _mp(material_properties)
+ {
+ }
+
+ bool computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ KelvinVector& sigma,
+ KelvinMatrix& C,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables) override;
+
+private:
+ MaterialProperties _mp;
+};
+
+} // namespace Solids
+} // namespace MaterialLib
+#include "Ehlers-impl.h"
+
+#endif // MATERIALLIB_SOLIDMODELS_EHLERS_H_
diff --git a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
index c2a74a9cb5ad11cc651f147fa9243f578123e622..4f8a1b4633d4fd31c1082eb35b4d1641c0dd07d4 100644
--- a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -13,6 +13,7 @@
#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
#include "MaterialLib/SolidModels/CreateLubby2.h"
+#include "MaterialLib/SolidModels/CreateEhlers.h"
#include "ProcessLib/Utils/ParseSecondaryVariables.h"
#include "SmallDeformationProcess.h"
@@ -74,7 +75,12 @@ createSmallDeformationProcess(
std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
material = nullptr;
- if (type == "LinearElasticIsotropic")
+ if (type == "Ehlers")
+ {
+ material = MaterialLib::Solids::createEhlers<DisplacementDim>(
+ parameters, constitutive_relation_config);
+ }
+ else if (type == "LinearElasticIsotropic")
{
material =
MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 18d3593d1ca889532428210f4cd18254746539dd..ace169d9adf2e6efaccf5e38e299cb0150cffb7f 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -112,6 +112,24 @@ struct SmallDeformationLocalAssemblerInterface
virtual std::vector<double> const& getIntPtSigmaYZ(
std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXX(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonZZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXY(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonXZ(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtEpsilonYZ(
+ std::vector<double>& cache) const = 0;
};
template <typename ShapeFunction, typename IntegrationMethod,
@@ -318,6 +336,44 @@ public:
return getIntPtSigma(cache, 5);
}
+ std::vector<double> const& getIntPtEpsilonXX(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 0);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 1);
+ }
+
+ std::vector<double> const& getIntPtEpsilonZZ(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 2);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXY(
+ std::vector<double>& cache) const override
+ {
+ return getIntPtEpsilon(cache, 3);
+ }
+
+ std::vector<double> const& getIntPtEpsilonXZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 4);
+ }
+
+ std::vector<double> const& getIntPtEpsilonYZ(
+ std::vector<double>& cache) const override
+ {
+ assert(DisplacementDim == 3);
+ return getIntPtEpsilon(cache, 5);
+ }
+
private:
std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
std::size_t const component) const
@@ -335,6 +391,20 @@ private:
return cache;
}
+ std::vector<double> const& getIntPtEpsilon(
+ std::vector<double>& cache, std::size_t const component) const
+ {
+ cache.clear();
+ cache.reserve(_ip_data.size());
+
+ for (auto const& ip_data : _ip_data) {
+ cache.push_back(ip_data._eps[component]);
+ }
+
+ return cache;
+ }
+
+
SmallDeformationProcessData<DisplacementDim>& _process_data;
std::vector<IntegrationPointData<BMatricesType, DisplacementDim>> _ip_data;
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index d21af9415860c629c5da91af1cc9c9ca85d4224d..59ebb0b628eded912339440a5d02e79aaa011b3a 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -116,6 +116,30 @@ private:
getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYZ));
}
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xx", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXX));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_yy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYY));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_zz", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonZZ));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "epsilon_xy", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXY));
}
void assembleConcreteProcess(const double t, GlobalVector const& x,
diff --git a/Tests/Data b/Tests/Data
index 6ffd7ec67ebd30d3c83754cf3127fd885e255528..74db8c483fcdd3231b40215a78a9357865a80f5e 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 6ffd7ec67ebd30d3c83754cf3127fd885e255528
+Subproject commit 74db8c483fcdd3231b40215a78a9357865a80f5e