diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/i_nonlinear_solver.md b/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/i_nonlinear_solver.md
deleted file mode 100644
index f6d09b0017c9833e58880901daacfa3631a13458..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/i_nonlinear_solver.md
+++ /dev/null
@@ -1,2 +0,0 @@
-Local nonlinear solver configuration used in the constitutive relation
-computation. The nonlinear solver is the NumLib::NewtonRaphson solver.
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_error_tolerance.md b/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_error_tolerance.md
deleted file mode 100644
index ccae39ea3ad802a7b785d4e758ad7b39757862cd..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_error_tolerance.md
+++ /dev/null
@@ -1 +0,0 @@
-\copydoc NumLib::NewtonRaphson::_tolerance_squared
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_maximum_iterations.md b/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_maximum_iterations.md
deleted file mode 100644
index c9b32edcfe5664c6a8bf684a2be2043ac4658c28..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/material/solid/constitutive_relation/Lubby2/nonlinear_solver/t_maximum_iterations.md
+++ /dev/null
@@ -1 +0,0 @@
-\copydoc NumLib::NewtonRaphson::_maximum_iterations
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/i_nonlinear_solver.md b/Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/i_nonlinear_solver.md
similarity index 100%
rename from Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/i_nonlinear_solver.md
rename to Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/i_nonlinear_solver.md
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/t_error_tolerance.md b/Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/t_error_tolerance.md
similarity index 100%
rename from Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/t_error_tolerance.md
rename to Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/t_error_tolerance.md
diff --git a/Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/t_maximum_iterations.md b/Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/t_maximum_iterations.md
similarity index 100%
rename from Documentation/ProjectFile/material/solid/constitutive_relation/Ehlers/nonlinear_solver/t_maximum_iterations.md
rename to Documentation/ProjectFile/material/solid/constitutive_relation/nonlinear_solver/t_maximum_iterations.md
diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..29de544b866b0f4f7a7628ec09786a25f616c97f
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
@@ -0,0 +1,75 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateConstitutiveRelation.cpp
+ * Created on July 10, 2018, 12:09 PM
+ */
+
+#include "CreateConstitutiveRelation.h"
+#include "CreateEhlers.h"
+#include "CreateLinearElasticIsotropic.h"
+#include "CreateLubby2.h"
+
+#include "MechanicsBase.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+createConstitutiveRelation(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ auto const constitutive_relation_config =
+ //! \ogs_file_param{material__solid__constitutive_relation}
+ config.getConfigSubtree("constitutive_relation");
+
+ auto const type =
+ //! \ogs_file_param{material__solid__constitutive_relation__type}
+ constitutive_relation_config.peekConfigParameter<std::string>("type");
+
+ if (type == "Ehlers")
+ {
+ return MaterialLib::Solids::Ehlers::createEhlers<DisplacementDim>(
+ parameters, constitutive_relation_config);
+ }
+ else if (type == "LinearElasticIsotropic")
+ {
+ return MaterialLib::Solids::createLinearElasticIsotropic<
+ DisplacementDim>(parameters, constitutive_relation_config);
+ }
+ else if (type == "Lubby2")
+ {
+ return MaterialLib::Solids::Lubby2::createLubby2<DisplacementDim>(
+ parameters, constitutive_relation_config);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Cannot construct constitutive relation of given type \'%s\'.",
+ type.c_str());
+ }
+}
+
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelation<2>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelation<3>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+}
+}
diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.h b/MaterialLib/SolidModels/CreateConstitutiveRelation.h
new file mode 100644
index 0000000000000000000000000000000000000000..15c43ebf81a20513d59894b0c8be6b0a036049b6
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateConstitutiveRelation.h
@@ -0,0 +1,50 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateConstitutiveRelation.h
+ * Created on July 10, 2018, 12:09 PM
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ProcessLib
+{
+struct ParameterBase;
+}
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+
+template <int DisplacementDim>
+std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+createConstitutiveRelation(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelation<2>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelation<3>(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+}
+}
diff --git a/MaterialLib/SolidModels/CreateEhlers.h b/MaterialLib/SolidModels/CreateEhlers.h
index 3969fef891a576da7202fbee04e4e01d34318ed0..7d181c1d85ec1439a4a23fdd277bfb80b959b8fa 100644
--- a/MaterialLib/SolidModels/CreateEhlers.h
+++ b/MaterialLib/SolidModels/CreateEhlers.h
@@ -9,6 +9,7 @@
#pragma once
+#include "CreateNewtonRaphsonSolverParameters.h"
#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
#include "Ehlers.h"
@@ -18,25 +19,6 @@ namespace Solids
{
namespace Ehlers
{
-inline NumLib::NewtonRaphsonSolverParameters
-createNewtonRaphsonSolverParameters(BaseLib::ConfigTree const& config)
-{
- DBUG("Create local nonlinear solver parameters.");
- auto const maximum_iterations =
- //! \ogs_file_param{material__solid__constitutive_relation__Ehlers__nonlinear_solver__maximum_iterations}
- config.getConfigParameter<int>("maximum_iterations");
-
- DBUG("\tmaximum_iterations: %d.", maximum_iterations);
-
- auto const error_tolerance =
- //! \ogs_file_param{material__solid__constitutive_relation__Ehlers__nonlinear_solver__error_tolerance}
- config.getConfigParameter<double>("error_tolerance");
-
- DBUG("\terror_tolerance: %g.", error_tolerance);
-
- return {maximum_iterations, error_tolerance};
-}
-
inline std::unique_ptr<DamagePropertiesParameters> createDamageProperties(
std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
BaseLib::ConfigTree const& config)
@@ -184,11 +166,8 @@ std::unique_ptr<SolidEhlers<DisplacementDim>> createEhlers(
createDamageProperties(parameters, *ehlers_damage_config);
}
- auto const& nonlinear_solver_config =
- //! \ogs_file_param{material__solid__constitutive_relation__Ehlers__nonlinear_solver}
- config.getConfigSubtree("nonlinear_solver");
auto const nonlinear_solver_parameters =
- createNewtonRaphsonSolverParameters(nonlinear_solver_config);
+ createNewtonRaphsonSolverParameters(config);
return std::make_unique<SolidEhlers<DisplacementDim>>(
nonlinear_solver_parameters, mp, std::move(ehlers_damage_properties));
diff --git a/MaterialLib/SolidModels/CreateLubby2.h b/MaterialLib/SolidModels/CreateLubby2.h
index 145d17e970a1ef847b09b299de0af8c15827bc1f..569d546e3805977d0838c53cf9a4aee6d51165c0 100644
--- a/MaterialLib/SolidModels/CreateLubby2.h
+++ b/MaterialLib/SolidModels/CreateLubby2.h
@@ -9,6 +9,7 @@
#pragma once
+#include "CreateNewtonRaphsonSolverParameters.h"
#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
#include "Lubby2.h"
@@ -19,25 +20,6 @@ namespace Solids
{
namespace Lubby2
{
-inline NumLib::NewtonRaphsonSolverParameters
-createNewtonRaphsonSolverParameters(BaseLib::ConfigTree const& config)
-{
- DBUG("Create local nonlinear solver parameters.");
- auto const maximum_iterations =
- //! \ogs_file_param{material__solid__constitutive_relation__Lubby2__nonlinear_solver__maximum_iterations}
- config.getConfigParameter<int>("maximum_iterations");
-
- DBUG("\tmaximum_iterations: %d.", maximum_iterations);
-
- auto const error_tolerance =
- //! \ogs_file_param{material__solid__constitutive_relation__Lubby2__nonlinear_solver__error_tolerance}
- config.getConfigParameter<double>("error_tolerance");
-
- DBUG("\terror_tolerance: %g.", error_tolerance);
-
- return {maximum_iterations, error_tolerance};
-}
-
template <int DisplacementDim>
std::unique_ptr<Lubby2<DisplacementDim>> createLubby2(
std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
@@ -117,11 +99,8 @@ std::unique_ptr<Lubby2<DisplacementDim>> createLubby2(
maxwell_viscosity, dependency_parameter_mK,
dependency_parameter_mvK, dependency_parameter_mvM};
- auto const& nonlinear_solver_config =
- //! \ogs_file_param{material__solid__constitutive_relation__Lubby2__nonlinear_solver}
- config.getConfigSubtree("nonlinear_solver");
auto const nonlinear_solver_parameters =
- createNewtonRaphsonSolverParameters(nonlinear_solver_config);
+ createNewtonRaphsonSolverParameters(config);
return std::unique_ptr<Lubby2<DisplacementDim>>{
new Lubby2<DisplacementDim>{nonlinear_solver_parameters, mp}};
diff --git a/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2ffd162f6fa36abce71c2bdcc9303fa04bc205d4
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp
@@ -0,0 +1,43 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateNewtonRaphsonSolverParameters.cpp
+ * Created on July 10, 2018, 11:32 AM
+ */
+
+#include "CreateNewtonRaphsonSolverParameters.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "NumLib/NewtonRaphson.h"
+
+namespace MaterialLib
+{
+NumLib::NewtonRaphsonSolverParameters createNewtonRaphsonSolverParameters(
+ BaseLib::ConfigTree const& config)
+{
+ DBUG("Create local nonlinear solver parameters.");
+ auto const& nonlinear_solver_config =
+ //! \ogs_file_param{material__solid__constitutive_relation__nonlinear_solver}
+ config.getConfigSubtree("nonlinear_solver");
+
+ auto const maximum_iterations =
+ //! \ogs_file_param{material__solid__constitutive_relation__nonlinear_solver__maximum_iterations}
+ nonlinear_solver_config.getConfigParameter<int>("maximum_iterations");
+
+ DBUG("\tmaximum_iterations: %d.", maximum_iterations);
+
+ auto const error_tolerance =
+ //! \ogs_file_param{material__solid__constitutive_relation__nonlinear_solver__error_tolerance}
+ nonlinear_solver_config.getConfigParameter<double>("error_tolerance");
+
+ DBUG("\terror_tolerance: %g.", error_tolerance);
+
+ return {maximum_iterations, error_tolerance};
+}
+}
diff --git a/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.h b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ecfb79d875d63b23b678d22f22f64590eb0f8da
--- /dev/null
+++ b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateNewtonRaphsonSolverParameters.h
+ * Created on July 10, 2018, 11:32 AM
+ */
+
+#pragma once
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace NumLib
+{
+struct NewtonRaphsonSolverParameters;
+}
+
+namespace MaterialLib
+{
+NumLib::NewtonRaphsonSolverParameters createNewtonRaphsonSolverParameters(
+ BaseLib::ConfigTree const& config);
+}
diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
deleted file mode 100644
index e6327db8434afaa341e42cc7d5db06f995c6c5aa..0000000000000000000000000000000000000000
--- a/MaterialLib/SolidModels/Ehlers-impl.h
+++ /dev/null
@@ -1,717 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-/**
- * 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
- *
- */
-#pragma once
-
-#include <boost/math/special_functions/pow.hpp>
-
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-namespace Ehlers
-{
-/// Special product of \c v with itself: \f$v \odot v\f$.
-/// The tensor \c v is given in Kelvin mapping.
-/// \note Implementation only for 2 and 3 dimensions.
-/// \attention Pay attention to the sign of the result, which normally would be
-/// negative, but the returned value is not negated. This has to do with \f$
-/// d(A^{-1})/dA = -A^{-1} \odot A^{-1} \f$.
-template <int DisplacementDim>
-MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> sOdotS(
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& v);
-
-template <int DisplacementDim>
-struct PhysicalStressWithInvariants final
-{
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- using KelvinVector =
- MathLib::KelvinVector::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)}
- {
- }
-
- KelvinVector value;
- KelvinVector D;
- double I_1;
- double J_2;
- double J_3;
-
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
-};
-
-/// 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(MaterialProperties const& mp,
- PhysicalStressWithInvariants<DisplacementDim> const& s,
- double const k)
-{
- 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 + mp.gamma * s.J_3 / (s.J_2 * std::sqrt(s.J_2)),
- mp.m) +
- mp.alpha / 2. * I_1_squared +
- boost::math::pow<2>(mp.delta) *
- boost::math::pow<2>(I_1_squared)) +
- mp.beta * s.I_1 + mp.epsilon * I_1_squared - k;
-}
-
-template <int DisplacementDim>
-typename SolidEhlers<DisplacementDim>::ResidualVectorType
-calculatePlasticResidual(
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_D,
- double const eps_V,
- PhysicalStressWithInvariants<DisplacementDim> const& s,
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_p_D,
- MathLib::KelvinVector::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,
- double const k,
- MaterialProperties const& mp)
-{
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- using KelvinVector =
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
-
- auto const& P_dev = Invariants::deviatoric_projection;
- auto const& identity2 = Invariants::identity2;
-
- double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
-
- typename SolidEhlers<DisplacementDim>::ResidualVectorType residual;
- // calculate stress residual
- residual.template segment<KelvinVectorSize>(0).noalias() =
- s.value / mp.G - 2 * (eps_D - eps_p_D) -
- mp.K / mp.G * (eps_V - eps_p_V) * identity2;
-
- // deviatoric plastic strain
- KelvinVector const sigma_D_inverse_D =
- P_dev * MathLib::KelvinVector::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{mp.gamma_p, theta, mp.m_p};
- double const sqrtPhi = std::sqrt(
- s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
- boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
- KelvinVector const flow_D = plasticFlowDeviatoricPart(
- s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.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 flow_V = plasticFlowVolumetricPart<DisplacementDim>(
- s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.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(mp, s, k) / mp.G;
- return residual;
-}
-
-template <int DisplacementDim>
-typename SolidEhlers<DisplacementDim>::JacobianMatrix calculatePlasticJacobian(
- double const dt,
- PhysicalStressWithInvariants<DisplacementDim> const& s,
- double const lambda,
- MaterialProperties const& mp)
-{
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- using KelvinVector =
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
- using KelvinMatrix =
- MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
-
- auto const& P_dev = Invariants::deviatoric_projection;
- auto const& identity2 = Invariants::identity2;
-
- double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
- OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.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 = MathLib::KelvinVector::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 + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
- boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
- KelvinVector const flow_D = plasticFlowDeviatoricPart(
- s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
- KelvinVector const lambda_flow_D = lambda * flow_D;
-
- typename SolidEhlers<DisplacementDim>::JacobianMatrix jacobian =
- SolidEhlers<DisplacementDim>::JacobianMatrix::Zero();
-
- // 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() = mp.K / mp.G * identity2;
-
- // G_14 and G_15 are zero
-
- // G_21 -- derivative of deviatoric flow
-
- double const gm_p = mp.gamma_p * mp.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) +
- (mp.alpha_p * s.I_1 +
- 4 * boost::math::pow<2>(mp.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 * sOdotS<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 - mp.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)) *
- mp.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 * mp.G *
- (-(mp.alpha_p * s.I_1 +
- 4 * boost::math::pow<2>(mp.delta_p) *
- boost::math::pow<3>(s.I_1)) /
- (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
- (mp.alpha_p * identity2 +
- 12 * boost::math::pow<2>(mp.delta_p * s.I_1) * identity2) /
- (2 * sqrtPhi) +
- 2 * mp.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, mp.alpha_p, mp.beta_p, mp.delta_p, mp.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, mp.m);
- double const gm = mp.gamma * mp.m;
- // derivative of yield function w.r.t. sigma
- KelvinVector const dF_dsigma =
- mp.G * (one_gt_pow_m * (s.D + gm * M0) +
- (mp.alpha * s.I_1 +
- 4 * boost::math::pow<2>(mp.delta) *
- boost::math::pow<3>(s.I_1)) *
- identity2) /
- (2. * sqrtPhi) +
- mp.G * (mp.beta + 2 * mp.epsilon_p * s.I_1) * identity2;
-
- jacobian
- .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
- .noalias() = dF_dsigma.transpose() / mp.G;
- }
-
- // G_54
- jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
- -mp.kappa * mp.hardening_coefficient / mp.G;
-
- // G_52, G_53, G_55 are zero
- return jacobian;
-}
-
-/// Calculates the derivative of the residuals with respect to total
-/// strain. Implementation fully implicit only.
-template <int DisplacementDim>
-MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> calculateDResidualDEps(
- double const K, double const G)
-{
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
-
- auto const& P_dev = Invariants::deviatoric_projection;
- auto const& P_sph = Invariants::spherical_projection;
- auto const& I =
- MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>::Identity();
-
- return -2. * I * P_dev - 3. * K / G * I * P_sph;
-}
-
-inline double calculateIsotropicHardening(double const kappa,
- double const hardening_coefficient,
- double const eps_p_eff)
-{
- return kappa * (1. + eps_p_eff * hardening_coefficient);
-}
-
-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 =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using Invariants = MathLib::KelvinVector::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;
-}
-
-/// Split the agglomerated solution vector in separate items. The arrangement
-/// must be the same as in the newton() function.
-template <typename ResidualVector, typename KelvinVector>
-std::tuple<KelvinVector, PlasticStrain<KelvinVector>, double>
-splitSolutionVector(ResidualVector const& solution)
-{
- static auto const size = KelvinVector::SizeAtCompileTime;
- return std::forward_as_tuple(
- solution.template segment<size>(size * 0),
- PlasticStrain<KelvinVector>{solution.template segment<size>(size * 1),
- solution[size * 2], solution[size * 2 + 1]},
- solution[size * 2 + 2]);
-}
-
-template <int DisplacementDim>
-boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
- std::unique_ptr<typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables>,
- typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
-SolidEhlers<DisplacementDim>::integrateStress(
- double const t,
- ProcessLib::SpatialPosition const& x,
- double const dt,
- KelvinVector const& eps_prev,
- KelvinVector const& eps,
- KelvinVector const& sigma_prev,
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
- material_state_variables)
-{
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &material_state_variables) != nullptr);
-
- StateVariables<DisplacementDim> state =
- static_cast<StateVariables<DisplacementDim> const&>(
- material_state_variables);
- state.setInitialConditions();
-
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
-
- // 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;
-
- // do the evaluation once per function call.
- MaterialProperties const mp(t, x, _mp);
-
- KelvinVector sigma = predict_sigma<DisplacementDim>(mp.G, mp.K, sigma_prev,
- eps, eps_prev, eps_V);
-
- KelvinMatrix tangentStiffness;
-
- PhysicalStressWithInvariants<DisplacementDim> s{mp.G * sigma};
- // Quit early if sigma is zero (nothing to do) or if we are still in elastic
- // zone.
- if ((sigma.squaredNorm() == 0 ||
- yieldFunction(
- mp, s,
- calculateIsotropicHardening(mp.kappa, mp.hardening_coefficient,
- state.eps_p.eff)) < 0))
- {
- tangentStiffness.setZero();
- tangentStiffness.template topLeftCorner<3, 3>().setConstant(
- mp.K - 2. / 3 * mp.G);
- tangentStiffness.noalias() += 2 * mp.G * KelvinMatrix::Identity();
- }
- else
- {
- // Linear solver for the newton loop is required after the loop with the
- // same matrix. This saves one decomposition.
- Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
- JacobianResidualSize, Eigen::RowMajor>>
- linear_solver;
-
- {
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<
- DisplacementDim>::value;
- using KelvinVector =
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
- using ResidualVectorType =
- Eigen::Matrix<double, JacobianResidualSize, 1>;
- using JacobianMatrix =
- Eigen::Matrix<double, JacobianResidualSize,
- JacobianResidualSize, Eigen::RowMajor>;
-
- JacobianMatrix jacobian;
-
- // Agglomerated solution vector construction. It is later split
- // into individual parts by splitSolutionVector().
- ResidualVectorType solution;
- solution << sigma, state.eps_p.D, state.eps_p.V, state.eps_p.eff, 0;
-
- auto const update_residual = [&](ResidualVectorType& residual) {
-
- auto const& eps_p_D =
- solution.template segment<KelvinVectorSize>(
- KelvinVectorSize);
- KelvinVector const eps_p_D_dot =
- (eps_p_D - state.eps_p_prev.D) / dt;
-
- double const& eps_p_V = solution[KelvinVectorSize * 2];
- double const eps_p_V_dot =
- (eps_p_V - state.eps_p_prev.V) / dt;
-
- double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
- double const eps_p_eff_dot =
- (eps_p_eff - state.eps_p_prev.eff) / dt;
-
- double const k_hardening = calculateIsotropicHardening(
- mp.kappa, mp.hardening_coefficient,
- solution[KelvinVectorSize * 2 + 1]);
- residual = calculatePlasticResidual<DisplacementDim>(
- eps_D, eps_V, s,
- solution.template segment<KelvinVectorSize>(
- KelvinVectorSize),
- eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
- eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
- k_hardening, mp);
- };
-
- auto const update_jacobian = [&](JacobianMatrix& jacobian) {
- jacobian = calculatePlasticJacobian<DisplacementDim>(
- dt, s, solution[KelvinVectorSize * 2 + 2], mp);
- };
-
- auto const update_solution =
- [&](ResidualVectorType const& increment) {
- solution += increment;
- s = PhysicalStressWithInvariants<DisplacementDim>{
- mp.G * solution.template segment<KelvinVectorSize>(0)};
- };
-
- auto newton_solver = NumLib::NewtonRaphson<
- decltype(linear_solver), JacobianMatrix,
- decltype(update_jacobian), ResidualVectorType,
- decltype(update_residual), decltype(update_solution)>(
- linear_solver, update_jacobian, update_residual,
- update_solution, _nonlinear_solver_parameters);
-
- auto const success_iterations = newton_solver.solve(jacobian);
-
- if (!success_iterations)
- return {};
-
- // 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);
-
- std::tie(sigma, state.eps_p, std::ignore) =
- splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
- }
-
- // 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>(mp.K, mp.G);
-
- tangentStiffness =
- mp.G *
- linear_solver.solve(-dresidual_deps)
- .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
- }
-
- KelvinVector sigma_final = mp.G * sigma;
-
- return {std::make_tuple(
- sigma_final,
- std::unique_ptr<
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
- new StateVariables<DisplacementDim>{
- static_cast<StateVariables<DisplacementDim> const&>(state)}},
- tangentStiffness)};
-}
-
-template <int DisplacementDim>
-std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
-SolidEhlers<DisplacementDim>::getInternalVariables() const
-{
- return {
- {"damage.kappa_d", 1,
- [](typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables const& state,
- std::vector<double>& cache) -> std::vector<double> const& {
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &state) != nullptr);
- auto const& ehlers_state =
- static_cast<StateVariables<DisplacementDim> const&>(state);
-
- cache.resize(1);
- cache.front() = ehlers_state.damage.kappa_d();
- return cache;
- }},
- {"damage.value", 1,
- [](typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables const& state,
- std::vector<double>& cache) -> std::vector<double> const& {
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &state) != nullptr);
- auto const& ehlers_state =
- static_cast<StateVariables<DisplacementDim> const&>(state);
-
- cache.resize(1);
- cache.front() = ehlers_state.damage.value();
- return cache;
- }},
- {"eps_p.D", KelvinVector::RowsAtCompileTime,
- [](typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables const& state,
- std::vector<double>& cache) -> std::vector<double> const& {
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &state) != nullptr);
- auto const& ehlers_state =
- static_cast<StateVariables<DisplacementDim> const&>(state);
-
- cache.resize(KelvinVector::RowsAtCompileTime);
- MathLib::toVector<KelvinVector>(cache,
- KelvinVector::RowsAtCompileTime) =
- MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
- ehlers_state.eps_p.D);
-
- return cache;
- }},
- {"eps_p.V", 1,
- [](typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables const& state,
- std::vector<double>& cache) -> std::vector<double> const& {
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &state) != nullptr);
- auto const& ehlers_state =
- static_cast<StateVariables<DisplacementDim> const&>(state);
-
- cache.resize(1);
- cache.front() = ehlers_state.eps_p.V;
- return cache;
- }},
- {"eps_p.eff", 1,
- [](typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables const& state,
- std::vector<double>& cache) -> std::vector<double> const& {
- assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
- &state) != nullptr);
- auto const& ehlers_state =
- static_cast<StateVariables<DisplacementDim> const&>(state);
-
- cache.resize(1);
- cache.front() = ehlers_state.eps_p.eff;
- return cache;
- }}};
-}
-
-} // namespace Ehlers
-} // namespace Solids
-} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
index 12a63de586f90a8e05c658179775dd93cac7f44f..99820eaf96f4b663939b88beca59dcbaeb2bfbe6 100644
--- a/MaterialLib/SolidModels/Ehlers.cpp
+++ b/MaterialLib/SolidModels/Ehlers.cpp
@@ -8,7 +8,33 @@
*/
#include "Ehlers.h"
-#include "Ehlers-impl.h"
+#include <boost/math/special_functions/pow.hpp>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+/**
+ * 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
+ *
+ */
namespace MaterialLib
{
@@ -16,6 +42,677 @@ namespace Solids
{
namespace Ehlers
{
+
+/// Special product of \c v with itself: \f$v \odot v\f$.
+/// The tensor \c v is given in Kelvin mapping.
+/// \note Implementation only for 2 and 3 dimensions.
+/// \attention Pay attention to the sign of the result, which normally would be
+/// negative, but the returned value is not negated. This has to do with \f$
+/// d(A^{-1})/dA = -A^{-1} \odot A^{-1} \f$.
+template <int DisplacementDim>
+MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> sOdotS(
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& v);
+
+template <int DisplacementDim>
+struct PhysicalStressWithInvariants final
+{
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ using KelvinVector =
+ MathLib::KelvinVector::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)}
+ {
+ }
+
+ KelvinVector value;
+ KelvinVector D;
+ double I_1;
+ double J_2;
+ double J_3;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+/// 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(MaterialProperties const& mp,
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ double const k)
+{
+ 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 + mp.gamma * s.J_3 / (s.J_2 * std::sqrt(s.J_2)),
+ mp.m) +
+ mp.alpha / 2. * I_1_squared +
+ boost::math::pow<2>(mp.delta) *
+ boost::math::pow<2>(I_1_squared)) +
+ mp.beta * s.I_1 + mp.epsilon * I_1_squared - k;
+}
+
+template <int DisplacementDim>
+typename SolidEhlers<DisplacementDim>::ResidualVectorType
+calculatePlasticResidual(
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_D,
+ double const eps_V,
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps_p_D,
+ MathLib::KelvinVector::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,
+ double const k,
+ MaterialProperties const& mp)
+{
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ using KelvinVector =
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& identity2 = Invariants::identity2;
+
+ double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
+
+ typename SolidEhlers<DisplacementDim>::ResidualVectorType residual;
+ // calculate stress residual
+ residual.template segment<KelvinVectorSize>(0).noalias() =
+ s.value / mp.G - 2 * (eps_D - eps_p_D) -
+ mp.K / mp.G * (eps_V - eps_p_V) * identity2;
+
+ // deviatoric plastic strain
+ KelvinVector const sigma_D_inverse_D =
+ P_dev * MathLib::KelvinVector::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{mp.gamma_p, theta, mp.m_p};
+ double const sqrtPhi = std::sqrt(
+ s.J_2 * one_gt.pow_m_p + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
+ boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
+ KelvinVector const flow_D = plasticFlowDeviatoricPart(
+ s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.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 flow_V = plasticFlowVolumetricPart<DisplacementDim>(
+ s, sqrtPhi, mp.alpha_p, mp.beta_p, mp.delta_p, mp.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(mp, s, k) / mp.G;
+ return residual;
+}
+
+template <int DisplacementDim>
+typename SolidEhlers<DisplacementDim>::JacobianMatrix calculatePlasticJacobian(
+ double const dt,
+ PhysicalStressWithInvariants<DisplacementDim> const& s,
+ double const lambda,
+ MaterialProperties const& mp)
+{
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ using KelvinVector =
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+ using KelvinMatrix =
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& identity2 = Invariants::identity2;
+
+ double const theta = s.J_3 / (s.J_2 * std::sqrt(s.J_2));
+ OnePlusGamma_pTheta const one_gt{mp.gamma_p, theta, mp.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 = MathLib::KelvinVector::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 + mp.alpha_p / 2. * boost::math::pow<2>(s.I_1) +
+ boost::math::pow<2>(mp.delta_p) * boost::math::pow<4>(s.I_1));
+ KelvinVector const flow_D = plasticFlowDeviatoricPart(
+ s, one_gt, sqrtPhi, dtheta_dsigma, mp.gamma_p, mp.m_p);
+ KelvinVector const lambda_flow_D = lambda * flow_D;
+
+ typename SolidEhlers<DisplacementDim>::JacobianMatrix jacobian =
+ SolidEhlers<DisplacementDim>::JacobianMatrix::Zero();
+
+ // 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() = mp.K / mp.G * identity2;
+
+ // G_14 and G_15 are zero
+
+ // G_21 -- derivative of deviatoric flow
+
+ double const gm_p = mp.gamma_p * mp.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) +
+ (mp.alpha_p * s.I_1 +
+ 4 * boost::math::pow<2>(mp.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 * sOdotS<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 - mp.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)) *
+ mp.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 * mp.G *
+ (-(mp.alpha_p * s.I_1 +
+ 4 * boost::math::pow<2>(mp.delta_p) *
+ boost::math::pow<3>(s.I_1)) /
+ (4 * boost::math::pow<3>(sqrtPhi)) * dPhi_dsigma +
+ (mp.alpha_p * identity2 +
+ 12 * boost::math::pow<2>(mp.delta_p * s.I_1) * identity2) /
+ (2 * sqrtPhi) +
+ 2 * mp.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, mp.alpha_p, mp.beta_p, mp.delta_p, mp.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, mp.m);
+ double const gm = mp.gamma * mp.m;
+ // derivative of yield function w.r.t. sigma
+ KelvinVector const dF_dsigma =
+ mp.G * (one_gt_pow_m * (s.D + gm * M0) +
+ (mp.alpha * s.I_1 +
+ 4 * boost::math::pow<2>(mp.delta) *
+ boost::math::pow<3>(s.I_1)) *
+ identity2) /
+ (2. * sqrtPhi) +
+ mp.G * (mp.beta + 2 * mp.epsilon_p * s.I_1) * identity2;
+
+ jacobian
+ .template block<1, KelvinVectorSize>(2 * KelvinVectorSize + 2, 0)
+ .noalias() = dF_dsigma.transpose() / mp.G;
+ }
+
+ // G_54
+ jacobian(2 * KelvinVectorSize + 2, 2 * KelvinVectorSize + 1) =
+ -mp.kappa * mp.hardening_coefficient / mp.G;
+
+ // G_52, G_53, G_55 are zero
+ return jacobian;
+}
+
+/// Calculates the derivative of the residuals with respect to total
+/// strain. Implementation fully implicit only.
+template <int DisplacementDim>
+MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> calculateDResidualDEps(
+ double const K, double const G)
+{
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+ auto const& P_dev = Invariants::deviatoric_projection;
+ auto const& P_sph = Invariants::spherical_projection;
+ auto const& I =
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>::Identity();
+
+ return -2. * I * P_dev - 3. * K / G * I * P_sph;
+}
+
+inline double calculateIsotropicHardening(double const kappa,
+ double const hardening_coefficient,
+ double const eps_p_eff)
+{
+ return kappa * (1. + eps_p_eff * hardening_coefficient);
+}
+
+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 =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using Invariants = MathLib::KelvinVector::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;
+}
+
+/// Split the agglomerated solution vector in separate items. The arrangement
+/// must be the same as in the newton() function.
+template <typename ResidualVector, typename KelvinVector>
+std::tuple<KelvinVector, PlasticStrain<KelvinVector>, double>
+splitSolutionVector(ResidualVector const& solution)
+{
+ static auto const size = KelvinVector::SizeAtCompileTime;
+ return std::forward_as_tuple(
+ solution.template segment<size>(size * 0),
+ PlasticStrain<KelvinVector>{solution.template segment<size>(size * 1),
+ solution[size * 2], solution[size * 2 + 1]},
+ solution[size * 2 + 2]);
+}
+
+template <int DisplacementDim>
+boost::optional<std::tuple<typename SolidEhlers<DisplacementDim>::KelvinVector,
+ std::unique_ptr<typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>,
+ typename SolidEhlers<DisplacementDim>::KelvinMatrix>>
+SolidEhlers<DisplacementDim>::integrateStress(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables)
+{
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &material_state_variables) != nullptr);
+
+ StateVariables<DisplacementDim> state =
+ static_cast<StateVariables<DisplacementDim> const&>(
+ material_state_variables);
+ state.setInitialConditions();
+
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+ // 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;
+
+ // do the evaluation once per function call.
+ MaterialProperties const mp(t, x, _mp);
+
+ KelvinVector sigma = predict_sigma<DisplacementDim>(mp.G, mp.K, sigma_prev,
+ eps, eps_prev, eps_V);
+
+ KelvinMatrix tangentStiffness;
+
+ PhysicalStressWithInvariants<DisplacementDim> s{mp.G * sigma};
+ // Quit early if sigma is zero (nothing to do) or if we are still in elastic
+ // zone.
+ if ((sigma.squaredNorm() == 0 ||
+ yieldFunction(
+ mp, s,
+ calculateIsotropicHardening(mp.kappa, mp.hardening_coefficient,
+ state.eps_p.eff)) < 0))
+ {
+ tangentStiffness.setZero();
+ tangentStiffness.template topLeftCorner<3, 3>().setConstant(
+ mp.K - 2. / 3 * mp.G);
+ tangentStiffness.noalias() += 2 * mp.G * KelvinMatrix::Identity();
+ }
+ else
+ {
+ // Linear solver for the newton loop is required after the loop with the
+ // same matrix. This saves one decomposition.
+ Eigen::FullPivLU<Eigen::Matrix<double, JacobianResidualSize,
+ JacobianResidualSize, Eigen::RowMajor>>
+ linear_solver;
+
+ {
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ using KelvinVector =
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
+ using ResidualVectorType =
+ Eigen::Matrix<double, JacobianResidualSize, 1>;
+ using JacobianMatrix =
+ Eigen::Matrix<double, JacobianResidualSize,
+ JacobianResidualSize, Eigen::RowMajor>;
+
+ JacobianMatrix jacobian;
+
+ // Agglomerated solution vector construction. It is later split
+ // into individual parts by splitSolutionVector().
+ ResidualVectorType solution;
+ solution << sigma, state.eps_p.D, state.eps_p.V, state.eps_p.eff, 0;
+
+ auto const update_residual = [&](ResidualVectorType& residual) {
+
+ auto const& eps_p_D =
+ solution.template segment<KelvinVectorSize>(
+ KelvinVectorSize);
+ KelvinVector const eps_p_D_dot =
+ (eps_p_D - state.eps_p_prev.D) / dt;
+
+ double const& eps_p_V = solution[KelvinVectorSize * 2];
+ double const eps_p_V_dot =
+ (eps_p_V - state.eps_p_prev.V) / dt;
+
+ double const& eps_p_eff = solution[KelvinVectorSize * 2 + 1];
+ double const eps_p_eff_dot =
+ (eps_p_eff - state.eps_p_prev.eff) / dt;
+
+ double const k_hardening = calculateIsotropicHardening(
+ mp.kappa, mp.hardening_coefficient,
+ solution[KelvinVectorSize * 2 + 1]);
+ residual = calculatePlasticResidual<DisplacementDim>(
+ eps_D, eps_V, s,
+ solution.template segment<KelvinVectorSize>(
+ KelvinVectorSize),
+ eps_p_D_dot, solution[KelvinVectorSize * 2], eps_p_V_dot,
+ eps_p_eff_dot, solution[KelvinVectorSize * 2 + 2],
+ k_hardening, mp);
+ };
+
+ auto const update_jacobian = [&](JacobianMatrix& jacobian) {
+ jacobian = calculatePlasticJacobian<DisplacementDim>(
+ dt, s, solution[KelvinVectorSize * 2 + 2], mp);
+ };
+
+ auto const update_solution =
+ [&](ResidualVectorType const& increment) {
+ solution += increment;
+ s = PhysicalStressWithInvariants<DisplacementDim>{
+ mp.G * solution.template segment<KelvinVectorSize>(0)};
+ };
+
+ auto newton_solver = NumLib::NewtonRaphson<
+ decltype(linear_solver), JacobianMatrix,
+ decltype(update_jacobian), ResidualVectorType,
+ decltype(update_residual), decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual,
+ update_solution, _nonlinear_solver_parameters);
+
+ auto const success_iterations = newton_solver.solve(jacobian);
+
+ if (!success_iterations)
+ return {};
+
+ // 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);
+
+ std::tie(sigma, state.eps_p, std::ignore) =
+ splitSolutionVector<ResidualVectorType, KelvinVector>(solution);
+ }
+
+ // 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>(mp.K, mp.G);
+
+ tangentStiffness =
+ mp.G *
+ linear_solver.solve(-dresidual_deps)
+ .template block<KelvinVectorSize, KelvinVectorSize>(0, 0);
+ }
+
+ KelvinVector sigma_final = mp.G * sigma;
+
+ return {std::make_tuple(
+ sigma_final,
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+ new StateVariables<DisplacementDim>{
+ static_cast<StateVariables<DisplacementDim> const&>(state)}},
+ tangentStiffness)};
+}
+
+template <int DisplacementDim>
+std::vector<typename MechanicsBase<DisplacementDim>::InternalVariable>
+SolidEhlers<DisplacementDim>::getInternalVariables() const
+{
+ return {
+ {"damage.kappa_d", 1,
+ [](typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const& state,
+ std::vector<double>& cache) -> std::vector<double> const& {
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &state) != nullptr);
+ auto const& ehlers_state =
+ static_cast<StateVariables<DisplacementDim> const&>(state);
+
+ cache.resize(1);
+ cache.front() = ehlers_state.damage.kappa_d();
+ return cache;
+ }},
+ {"damage.value", 1,
+ [](typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const& state,
+ std::vector<double>& cache) -> std::vector<double> const& {
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &state) != nullptr);
+ auto const& ehlers_state =
+ static_cast<StateVariables<DisplacementDim> const&>(state);
+
+ cache.resize(1);
+ cache.front() = ehlers_state.damage.value();
+ return cache;
+ }},
+ {"eps_p.D", KelvinVector::RowsAtCompileTime,
+ [](typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const& state,
+ std::vector<double>& cache) -> std::vector<double> const& {
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &state) != nullptr);
+ auto const& ehlers_state =
+ static_cast<StateVariables<DisplacementDim> const&>(state);
+
+ cache.resize(KelvinVector::RowsAtCompileTime);
+ MathLib::toVector<KelvinVector>(cache,
+ KelvinVector::RowsAtCompileTime) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ ehlers_state.eps_p.D);
+
+ return cache;
+ }},
+ {"eps_p.V", 1,
+ [](typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const& state,
+ std::vector<double>& cache) -> std::vector<double> const& {
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &state) != nullptr);
+ auto const& ehlers_state =
+ static_cast<StateVariables<DisplacementDim> const&>(state);
+
+ cache.resize(1);
+ cache.front() = ehlers_state.eps_p.V;
+ return cache;
+ }},
+ {"eps_p.eff", 1,
+ [](typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables const& state,
+ std::vector<double>& cache) -> std::vector<double> const& {
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &state) != nullptr);
+ auto const& ehlers_state =
+ static_cast<StateVariables<DisplacementDim> const&>(state);
+
+ cache.resize(1);
+ cache.front() = ehlers_state.eps_p.eff;
+ return cache;
+ }}};
+}
+
template class SolidEhlers<2>;
template class SolidEhlers<3>;
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic-impl.h b/MaterialLib/SolidModels/LinearElasticIsotropic-impl.h
deleted file mode 100644
index f7b0e7ea3e3d37d8d80e9071eebd8de2c6bae614..0000000000000000000000000000000000000000
--- a/MaterialLib/SolidModels/LinearElasticIsotropic-impl.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "LinearElasticIsotropic.h"
-
-namespace MaterialLib
-{
-namespace Solids
-{
-template <int DisplacementDim>
-boost::optional<
- std::tuple<typename LinearElasticIsotropic<DisplacementDim>::KelvinVector,
- std::unique_ptr<typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables>,
- typename LinearElasticIsotropic<DisplacementDim>::KelvinMatrix>>
-LinearElasticIsotropic<DisplacementDim>::integrateStress(
- double const t,
- ProcessLib::SpatialPosition const& x,
- double const /*dt*/,
- KelvinVector const& eps_prev,
- KelvinVector const& eps,
- KelvinVector const& sigma_prev,
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
- material_state_variables)
-{
- KelvinMatrix C = KelvinMatrix::Zero();
-
- C.template topLeftCorner<3, 3>().setConstant(_mp.lambda(t, x));
- C.noalias() += 2 * _mp.mu(t, x) * KelvinMatrix::Identity();
-
- KelvinVector sigma = sigma_prev + C * (eps - eps_prev);
-
- return {std::make_tuple(
- sigma,
- std::unique_ptr<
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
- new MaterialStateVariables{
- static_cast<MaterialStateVariables const&>(
- material_state_variables)}},
- C)};
-}
-
-extern template class LinearElasticIsotropic<2>;
-extern template class LinearElasticIsotropic<3>;
-
-} // namespace Solids
-} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
index 1e5d988b4123a999e969dbfad2fda98927aade03..9aa837a7f7b711db90bf9efbe03cd54fdab24bc3 100644
--- a/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
+++ b/MaterialLib/SolidModels/LinearElasticIsotropic.cpp
@@ -8,12 +8,43 @@
*/
#include "LinearElasticIsotropic.h"
-#include "LinearElasticIsotropic-impl.h"
namespace MaterialLib
{
namespace Solids
{
+template <int DisplacementDim>
+boost::optional<
+ std::tuple<typename LinearElasticIsotropic<DisplacementDim>::KelvinVector,
+ std::unique_ptr<typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>,
+ typename LinearElasticIsotropic<DisplacementDim>::KelvinMatrix>>
+LinearElasticIsotropic<DisplacementDim>::integrateStress(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const /*dt*/,
+ KelvinVector const& eps_prev,
+ KelvinVector const& eps,
+ KelvinVector const& sigma_prev,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables)
+{
+ KelvinMatrix C = KelvinMatrix::Zero();
+
+ C.template topLeftCorner<3, 3>().setConstant(_mp.lambda(t, x));
+ C.noalias() += 2 * _mp.mu(t, x) * KelvinMatrix::Identity();
+
+ KelvinVector sigma = sigma_prev + C * (eps - eps_prev);
+
+ return {std::make_tuple(
+ sigma,
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+ new MaterialStateVariables{
+ static_cast<MaterialStateVariables const&>(
+ material_state_variables)}},
+ C)};
+}
template class LinearElasticIsotropic<2>;
template class LinearElasticIsotropic<3>;
diff --git a/MaterialLib/SolidModels/Lubby2-impl.h b/MaterialLib/SolidModels/Lubby2-impl.h
deleted file mode 100644
index 6a3ce0fbe70c903c3f36dd02af5a27cf102378b2..0000000000000000000000000000000000000000
--- a/MaterialLib/SolidModels/Lubby2-impl.h
+++ /dev/null
@@ -1,306 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-namespace MaterialLib
-{
-namespace Solids
-{
-namespace Lubby2
-{
-/// Calculates the 18x6 derivative of the residuals with respect to total
-/// strain.
-///
-/// Function definition can not be moved into implementation because of a
-/// MSVC compiler errors. See
-/// http://stackoverflow.com/questions/1484885/strange-vc-compile-error-c2244
-/// and https://support.microsoft.com/en-us/kb/930198
-template <int DisplacementDim>
-Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
- Lubby2<DisplacementDim>::KelvinVectorSize>
-calculatedGdEBurgers()
-{
- Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
- Lubby2<DisplacementDim>::KelvinVectorSize>
- dGdE =
- Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
- Lubby2<DisplacementDim>::KelvinVectorSize>::Zero();
- dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
- Lubby2<DisplacementDim>::KelvinVectorSize>()
- .diagonal()
- .setConstant(-2);
- return dGdE;
-}
-
-template <int DisplacementDim, typename LinearSolver>
-MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> tangentStiffnessA(
- double const GM0, double const KM0, LinearSolver const& linear_solver)
-{
- // Calculate dGdE for time step
- auto const dGdE = calculatedGdEBurgers<DisplacementDim>();
-
- // Consistent tangent from local Newton iteration of material
- // functionals.
- // Only the upper left block is relevant for the global tangent.
- static int const KelvinVectorSize =
- MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
- using KelvinMatrix =
- MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
-
- KelvinMatrix const dzdE =
- linear_solver.solve(-dGdE)
- .template topLeftCorner<KelvinVectorSize, KelvinVectorSize>();
-
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- auto const& P_sph = Invariants::spherical_projection;
- auto const& P_dev = Invariants::deviatoric_projection;
-
- KelvinMatrix C = GM0 * dzdE * P_dev + 3. * KM0 * P_sph;
- return C;
-};
-
-template <int DisplacementDim>
-boost::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
- std::unique_ptr<typename MechanicsBase<
- DisplacementDim>::MaterialStateVariables>,
- typename Lubby2<DisplacementDim>::KelvinMatrix>>
-Lubby2<DisplacementDim>::integrateStress(
- double const t,
- ProcessLib::SpatialPosition const& x,
- double const dt,
- KelvinVector const& /*eps_prev*/,
- KelvinVector const& eps,
- KelvinVector const& /*sigma_prev*/,
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
- material_state_variables)
-{
- using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
-
- assert(dynamic_cast<MaterialStateVariables const*>(
- &material_state_variables) != nullptr);
- MaterialStateVariables state(
- static_cast<MaterialStateVariables const&>(material_state_variables));
- state.setInitialConditions();
-
- auto local_lubby2_properties =
- detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
-
- // calculation of deviatoric parts
- auto const& P_dev = Invariants::deviatoric_projection;
- KelvinVector const epsd_i = P_dev * eps;
-
- // initial guess as elastic predictor
- KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
-
- // Calculate effective stress and update material properties
- double sig_eff = Invariants::equivalentStress(sigd_j);
- local_lubby2_properties.update(sig_eff);
-
- using LocalJacobianMatrix =
- Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
- Eigen::RowMajor>;
-
- // Linear solver for the newton loop is required after the loop with the
- // same matrix. This saves one decomposition.
- Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
-
- // Different solvers are available for the solution of the local system.
- // TODO Make the following choice of linear solvers available from the
- // input file configuration:
- // K_loc.partialPivLu().solve(-res_loc);
- // K_loc.fullPivLu().solve(-res_loc);
- // K_loc.householderQr().solve(-res_loc);
- // K_loc.colPivHouseholderQr().solve(res_loc);
- // K_loc.fullPivHouseholderQr().solve(-res_loc);
- // K_loc.llt().solve(-res_loc);
- // K_loc.ldlt().solve(-res_loc);
-
- LocalJacobianMatrix K_loc;
- { // Local Newton solver
- using LocalResidualVector =
- Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
-
- auto const update_residual = [&](LocalResidualVector& residual) {
- calculateResidualBurgers(
- dt, epsd_i, sigd_j, state.eps_K_j, state.eps_K_t, state.eps_M_j,
- state.eps_M_t, residual, local_lubby2_properties);
- };
-
- auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
- calculateJacobianBurgers(
- t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
- local_lubby2_properties); // for solution dependent Jacobians
- };
-
- auto const update_solution = [&](LocalResidualVector const& increment) {
- // increment solution vectors
- sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
- KelvinVectorSize * 0);
- state.eps_K_j.noalias() +=
- increment.template segment<KelvinVectorSize>(KelvinVectorSize *
- 1);
- state.eps_M_j.noalias() +=
- increment.template segment<KelvinVectorSize>(KelvinVectorSize *
- 2);
-
- // Calculate effective stress and update material properties
- sig_eff = MathLib::KelvinVector::Invariants<
- KelvinVectorSize>::equivalentStress(sigd_j);
- local_lubby2_properties.update(sig_eff);
- };
-
- auto newton_solver = NumLib::NewtonRaphson<
- decltype(linear_solver), LocalJacobianMatrix,
- decltype(update_jacobian), LocalResidualVector,
- decltype(update_residual), decltype(update_solution)>(
- linear_solver, update_jacobian, update_residual, update_solution,
- _nonlinear_solver_parameters);
-
- auto const success_iterations = newton_solver.solve(K_loc);
-
- if (!success_iterations)
- return {};
-
- // 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(K_loc);
- }
-
- KelvinMatrix C =
- tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
- local_lubby2_properties.KM0,
- linear_solver);
-
- // Hydrostatic part for the stress and the tangent.
- double const eps_i_trace = Invariants::trace(eps);
- KelvinVector const sigma =
- local_lubby2_properties.GM0 * sigd_j +
- local_lubby2_properties.KM0 * eps_i_trace * Invariants::identity2;
- return {std::make_tuple(
- sigma,
- std::unique_ptr<
- typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
- new MaterialStateVariables{state}},
- C)};
-}
-
-template <int DisplacementDim>
-void Lubby2<DisplacementDim>::calculateResidualBurgers(
- const double dt,
- const KelvinVector& strain_curr,
- const KelvinVector& stress_curr,
- KelvinVector& strain_Kel_curr,
- const KelvinVector& strain_Kel_t,
- KelvinVector& strain_Max_curr,
- const KelvinVector& strain_Max_t,
- ResidualVector& res,
- detail::LocalLubby2Properties<DisplacementDim> const& properties)
-{
- // calculate stress residual
- res.template segment<KelvinVectorSize>(0).noalias() =
- stress_curr - 2. * (strain_curr - strain_Kel_curr - strain_Max_curr);
-
- // calculate Kelvin strain residual
- res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
- 1. / dt * (strain_Kel_curr - strain_Kel_t) -
- 1. / (2. * properties.etaK) * (properties.GM0 * stress_curr -
- 2. * properties.GK * strain_Kel_curr);
-
- // calculate Maxwell strain residual
- res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
- 1. / dt * (strain_Max_curr - strain_Max_t) -
- 0.5 * properties.GM0 / properties.etaM * stress_curr;
-}
-
-template <int DisplacementDim>
-void Lubby2<DisplacementDim>::calculateJacobianBurgers(
- double const t,
- ProcessLib::SpatialPosition const& x,
- const double dt,
- JacobianMatrix& Jac,
- double s_eff,
- const KelvinVector& sig_i,
- const KelvinVector& eps_K_i,
- detail::LocalLubby2Properties<DisplacementDim> const& properties)
-{
- Jac.setZero();
-
- // build G_11
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0).setIdentity();
-
- // build G_12
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
- .diagonal()
- .setConstant(2);
-
- // build G_13
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
- 2 * KelvinVectorSize)
- .diagonal()
- .setConstant(2);
-
- // build G_21
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
- .noalias() =
- -0.5 * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
- if (s_eff > 0.)
- {
- KelvinVector const eps_K_aid =
- 1. / (properties.etaK * properties.etaK) *
- (properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
-
- KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
- properties.GM0 / s_eff * sig_i;
- KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
- properties.etaK / s_eff * sig_i;
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
- 0)
- .noalias() += 0.5 * eps_K_aid * dmu_vK.transpose() +
- 1. / properties.etaK * eps_K_i * dG_K.transpose();
- }
-
- // build G_22
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
- KelvinVectorSize)
- .diagonal()
- .setConstant(1. / dt + properties.GK / properties.etaK);
-
- // nothing to do for G_23
-
- // build G_31
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
- 0)
- .noalias() =
- -0.5 * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
- if (s_eff > 0.)
- {
- KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
- properties.etaM / s_eff * sig_i;
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(
- 2 * KelvinVectorSize, 0)
- .noalias() += 0.5 * properties.GM0 /
- (properties.etaM * properties.etaM) * sig_i *
- dmu_vM.transpose();
- }
-
- // nothing to do for G_32
-
- // build G_33
- Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
- 2 * KelvinVectorSize)
- .diagonal()
- .setConstant(1. / dt);
-}
-
-} // namespace Lubby2
-} // namespace Solids
-} // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/Lubby2.cpp b/MaterialLib/SolidModels/Lubby2.cpp
index 983ab6e48d05cc00ce550a3053bc30cb8306970f..97f40ac3bca4550a5791a7cb20257914197c9eeb 100644
--- a/MaterialLib/SolidModels/Lubby2.cpp
+++ b/MaterialLib/SolidModels/Lubby2.cpp
@@ -8,7 +8,6 @@
*/
#include "Lubby2.h"
-#include "Lubby2-impl.h"
namespace MaterialLib
{
@@ -16,6 +15,293 @@ namespace Solids
{
namespace Lubby2
{
+
+/// Calculates the 18x6 derivative of the residuals with respect to total
+/// strain.
+///
+/// Function definition can not be moved into implementation because of a
+/// MSVC compiler errors. See
+/// http://stackoverflow.com/questions/1484885/strange-vc-compile-error-c2244
+/// and https://support.microsoft.com/en-us/kb/930198
+template <int DisplacementDim>
+Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
+ Lubby2<DisplacementDim>::KelvinVectorSize>
+calculatedGdEBurgers()
+{
+ Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
+ Lubby2<DisplacementDim>::KelvinVectorSize>
+ dGdE =
+ Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
+ Lubby2<DisplacementDim>::KelvinVectorSize>::Zero();
+ dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
+ Lubby2<DisplacementDim>::KelvinVectorSize>()
+ .diagonal()
+ .setConstant(-2);
+ return dGdE;
+}
+
+template <int DisplacementDim, typename LinearSolver>
+MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> tangentStiffnessA(
+ double const GM0, double const KM0, LinearSolver const& linear_solver)
+{
+ // Calculate dGdE for time step
+ auto const dGdE = calculatedGdEBurgers<DisplacementDim>();
+
+ // Consistent tangent from local Newton iteration of material
+ // functionals.
+ // Only the upper left block is relevant for the global tangent.
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
+ using KelvinMatrix =
+ MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
+
+ KelvinMatrix const dzdE =
+ linear_solver.solve(-dGdE)
+ .template topLeftCorner<KelvinVectorSize, KelvinVectorSize>();
+
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ auto const& P_sph = Invariants::spherical_projection;
+ auto const& P_dev = Invariants::deviatoric_projection;
+
+ KelvinMatrix C = GM0 * dzdE * P_dev + 3. * KM0 * P_sph;
+ return C;
+};
+
+template <int DisplacementDim>
+boost::optional<std::tuple<typename Lubby2<DisplacementDim>::KelvinVector,
+ std::unique_ptr<typename MechanicsBase<
+ DisplacementDim>::MaterialStateVariables>,
+ typename Lubby2<DisplacementDim>::KelvinMatrix>>
+Lubby2<DisplacementDim>::integrateStress(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const dt,
+ KelvinVector const& /*eps_prev*/,
+ KelvinVector const& eps,
+ KelvinVector const& /*sigma_prev*/,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables const&
+ material_state_variables)
+{
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+
+ assert(dynamic_cast<MaterialStateVariables const*>(
+ &material_state_variables) != nullptr);
+ MaterialStateVariables state(
+ static_cast<MaterialStateVariables const&>(material_state_variables));
+ state.setInitialConditions();
+
+ auto local_lubby2_properties =
+ detail::LocalLubby2Properties<DisplacementDim>{t, x, _mp};
+
+ // calculation of deviatoric parts
+ auto const& P_dev = Invariants::deviatoric_projection;
+ KelvinVector const epsd_i = P_dev * eps;
+
+ // initial guess as elastic predictor
+ KelvinVector sigd_j = 2.0 * (epsd_i - state.eps_M_t - state.eps_K_t);
+
+ // Calculate effective stress and update material properties
+ double sig_eff = Invariants::equivalentStress(sigd_j);
+ local_lubby2_properties.update(sig_eff);
+
+ using LocalJacobianMatrix =
+ Eigen::Matrix<double, KelvinVectorSize * 3, KelvinVectorSize * 3,
+ Eigen::RowMajor>;
+
+ // Linear solver for the newton loop is required after the loop with the
+ // same matrix. This saves one decomposition.
+ Eigen::FullPivLU<LocalJacobianMatrix> linear_solver;
+
+ // Different solvers are available for the solution of the local system.
+ // TODO Make the following choice of linear solvers available from the
+ // input file configuration:
+ // K_loc.partialPivLu().solve(-res_loc);
+ // K_loc.fullPivLu().solve(-res_loc);
+ // K_loc.householderQr().solve(-res_loc);
+ // K_loc.colPivHouseholderQr().solve(res_loc);
+ // K_loc.fullPivHouseholderQr().solve(-res_loc);
+ // K_loc.llt().solve(-res_loc);
+ // K_loc.ldlt().solve(-res_loc);
+
+ LocalJacobianMatrix K_loc;
+ { // Local Newton solver
+ using LocalResidualVector =
+ Eigen::Matrix<double, KelvinVectorSize * 3, 1>;
+
+ auto const update_residual = [&](LocalResidualVector& residual) {
+ calculateResidualBurgers(
+ dt, epsd_i, sigd_j, state.eps_K_j, state.eps_K_t, state.eps_M_j,
+ state.eps_M_t, residual, local_lubby2_properties);
+ };
+
+ auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
+ calculateJacobianBurgers(
+ t, x, dt, jacobian, sig_eff, sigd_j, state.eps_K_j,
+ local_lubby2_properties); // for solution dependent Jacobians
+ };
+
+ auto const update_solution = [&](LocalResidualVector const& increment) {
+ // increment solution vectors
+ sigd_j.noalias() += increment.template segment<KelvinVectorSize>(
+ KelvinVectorSize * 0);
+ state.eps_K_j.noalias() +=
+ increment.template segment<KelvinVectorSize>(KelvinVectorSize *
+ 1);
+ state.eps_M_j.noalias() +=
+ increment.template segment<KelvinVectorSize>(KelvinVectorSize *
+ 2);
+
+ // Calculate effective stress and update material properties
+ sig_eff = MathLib::KelvinVector::Invariants<
+ KelvinVectorSize>::equivalentStress(sigd_j);
+ local_lubby2_properties.update(sig_eff);
+ };
+
+ auto newton_solver = NumLib::NewtonRaphson<
+ decltype(linear_solver), LocalJacobianMatrix,
+ decltype(update_jacobian), LocalResidualVector,
+ decltype(update_residual), decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual, update_solution,
+ _nonlinear_solver_parameters);
+
+ auto const success_iterations = newton_solver.solve(K_loc);
+
+ if (!success_iterations)
+ return {};
+
+ // 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(K_loc);
+ }
+
+ KelvinMatrix C =
+ tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
+ local_lubby2_properties.KM0,
+ linear_solver);
+
+ // Hydrostatic part for the stress and the tangent.
+ double const eps_i_trace = Invariants::trace(eps);
+ KelvinVector const sigma =
+ local_lubby2_properties.GM0 * sigd_j +
+ local_lubby2_properties.KM0 * eps_i_trace * Invariants::identity2;
+ return {std::make_tuple(
+ sigma,
+ std::unique_ptr<
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables>{
+ new MaterialStateVariables{state}},
+ C)};
+}
+
+template <int DisplacementDim>
+void Lubby2<DisplacementDim>::calculateResidualBurgers(
+ const double dt,
+ const KelvinVector& strain_curr,
+ const KelvinVector& stress_curr,
+ KelvinVector& strain_Kel_curr,
+ const KelvinVector& strain_Kel_t,
+ KelvinVector& strain_Max_curr,
+ const KelvinVector& strain_Max_t,
+ ResidualVector& res,
+ detail::LocalLubby2Properties<DisplacementDim> const& properties)
+{
+ // calculate stress residual
+ res.template segment<KelvinVectorSize>(0).noalias() =
+ stress_curr - 2. * (strain_curr - strain_Kel_curr - strain_Max_curr);
+
+ // calculate Kelvin strain residual
+ res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
+ 1. / dt * (strain_Kel_curr - strain_Kel_t) -
+ 1. / (2. * properties.etaK) * (properties.GM0 * stress_curr -
+ 2. * properties.GK * strain_Kel_curr);
+
+ // calculate Maxwell strain residual
+ res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
+ 1. / dt * (strain_Max_curr - strain_Max_t) -
+ 0.5 * properties.GM0 / properties.etaM * stress_curr;
+}
+
+template <int DisplacementDim>
+void Lubby2<DisplacementDim>::calculateJacobianBurgers(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ const double dt,
+ JacobianMatrix& Jac,
+ double s_eff,
+ const KelvinVector& sig_i,
+ const KelvinVector& eps_K_i,
+ detail::LocalLubby2Properties<DisplacementDim> const& properties)
+{
+ Jac.setZero();
+
+ // build G_11
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0).setIdentity();
+
+ // build G_12
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
+ .diagonal()
+ .setConstant(2);
+
+ // build G_13
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
+ 2 * KelvinVectorSize)
+ .diagonal()
+ .setConstant(2);
+
+ // build G_21
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
+ .noalias() =
+ -0.5 * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
+ if (s_eff > 0.)
+ {
+ KelvinVector const eps_K_aid =
+ 1. / (properties.etaK * properties.etaK) *
+ (properties.GM0 * sig_i - 2. * properties.GK * eps_K_i);
+
+ KelvinVector const dG_K = 1.5 * _mp.mK(t, x)[0] * properties.GK *
+ properties.GM0 / s_eff * sig_i;
+ KelvinVector const dmu_vK = 1.5 * _mp.mvK(t, x)[0] * properties.GM0 *
+ properties.etaK / s_eff * sig_i;
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
+ 0)
+ .noalias() += 0.5 * eps_K_aid * dmu_vK.transpose() +
+ 1. / properties.etaK * eps_K_i * dG_K.transpose();
+ }
+
+ // build G_22
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
+ KelvinVectorSize)
+ .diagonal()
+ .setConstant(1. / dt + properties.GK / properties.etaK);
+
+ // nothing to do for G_23
+
+ // build G_31
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
+ 0)
+ .noalias() =
+ -0.5 * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
+ if (s_eff > 0.)
+ {
+ KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
+ properties.etaM / s_eff * sig_i;
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(
+ 2 * KelvinVectorSize, 0)
+ .noalias() += 0.5 * properties.GM0 /
+ (properties.etaM * properties.etaM) * sig_i *
+ dmu_vM.transpose();
+ }
+
+ // nothing to do for G_32
+
+ // build G_33
+ Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
+ 2 * KelvinVectorSize)
+ .diagonal()
+ .setConstant(1. / dt);
+}
+
template class Lubby2<2>;
template class Lubby2<3>;
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 61782e91cc3a90717f19fa5da681e41dee02bc0a..4914d6fce61003187e7b618e43e4ed74b5bbd158 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -11,8 +11,10 @@
#include <cassert>
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
#include "HydroMechanicsProcess.h"
#include "HydroMechanicsProcessData.h"
@@ -100,29 +102,9 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
}
// Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
-
- auto const type =
- //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material = nullptr;
- if (type == "LinearElasticIsotropic")
- {
- material =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
+ auto material =
+ MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+ parameters, config);
// Intrinsic permeability
auto& intrinsic_permeability = findParameter<double>(
diff --git a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 44a21ff9b3e481b3a9d52c2fc173ce8e48d8c64e..5a8552fb5177c29132a0c66a5d615db28206bdbe 100644
--- a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -14,7 +14,7 @@
#include "MaterialLib/FractureModels/CreateCohesiveZoneModeI.h"
#include "MaterialLib/FractureModels/CreateLinearElasticIsotropic.h"
#include "MaterialLib/FractureModels/CreateMohrCoulomb.h"
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
@@ -125,28 +125,9 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
// Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
-
- auto const type =
- //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS_WITH_LIE__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<GlobalDim>> material =
- nullptr;
- if (type == "LinearElasticIsotropic")
- {
- material = MaterialLib::Solids::createLinearElasticIsotropic<GlobalDim>(
- parameters, constitutive_relation_config);
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
+ auto material =
+ MaterialLib::Solids::createConstitutiveRelation<GlobalDim>(
+ parameters, config);
// Intrinsic permeability
auto& intrinsic_permeability = findParameter<double>(
diff --git a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
index 9038979f191a1a1a36c09551d5c871c986930339..b834695aca4de61eb4f74431dd8c51fcdfb81463 100644
--- a/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -14,7 +14,7 @@
#include "MaterialLib/FractureModels/CreateCohesiveZoneModeI.h"
#include "MaterialLib/FractureModels/CreateLinearElasticIsotropic.h"
#include "MaterialLib/FractureModels/CreateMohrCoulomb.h"
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
@@ -99,29 +99,9 @@ std::unique_ptr<Process> createSmallDeformationProcess(
process_variables.push_back(std::move(per_process_variables));
// Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
-
- auto const type =
- //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION_WITH_LIE__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material = nullptr;
- if (type == "LinearElasticIsotropic")
- {
- material =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
+ auto material =
+ MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+ parameters, config);
// Fracture constitutive relation.
// read type;
diff --git a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
index 0df0542d83eeee8e96ef7f1a30631d2a8369ef91..ca1ee9024aba05f118b75ed5b3afb62f6ec2c533 100644
--- a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -11,10 +11,9 @@
#include <cassert>
-#include "MaterialLib/SolidModels/CreateEhlers.h"
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
-#include "MaterialLib/SolidModels/CreateLubby2.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
#include "SmallDeformationProcess.h"
#include "SmallDeformationProcessData.h"
@@ -66,39 +65,9 @@ createSmallDeformationProcess(
process_variables.push_back(std::move(per_process_variables));
// Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
-
- auto const type =
- //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material = nullptr;
- if (type == "Ehlers")
- {
- material = MaterialLib::Solids::Ehlers::createEhlers<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else if (type == "LinearElasticIsotropic")
- {
- material =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else if (type == "Lubby2")
- {
- material = MaterialLib::Solids::Lubby2::createLubby2<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
+ auto material =
+ MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+ parameters, config);
// Solid density
auto& solid_density = findParameter<double>(
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index e8994170c3b507c1e5c43fbfd4618d8b8bc5c89a..e57ed7bf8cc78562db7ea3a9c01a68fc6952bd08 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -11,10 +11,10 @@
#include <cassert>
-#include "MaterialLib/SolidModels/CreateEhlers.h"
-#include "MaterialLib/SolidModels/CreateLinearElasticIsotropic.h"
-#include "MaterialLib/SolidModels/CreateLubby2.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
#include "ThermoMechanicsProcess.h"
#include "ThermoMechanicsProcessData.h"
@@ -101,39 +101,9 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
}
// Constitutive relation.
- // read type;
- auto const constitutive_relation_config =
- //! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__constitutive_relation}
- config.getConfigSubtree("constitutive_relation");
-
- auto const type =
- //! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__constitutive_relation__type}
- constitutive_relation_config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
- material = nullptr;
- if (type == "Ehlers")
- {
- material = MaterialLib::Solids::Ehlers::createEhlers<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else if (type == "LinearElasticIsotropic")
- {
- material =
- MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else if (type == "Lubby2")
- {
- material = MaterialLib::Solids::Lubby2::createLubby2<DisplacementDim>(
- parameters, constitutive_relation_config);
- }
- else
- {
- OGS_FATAL(
- "Cannot construct constitutive relation of given type \'%s\'.",
- type.c_str());
- }
+ auto material =
+ MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+ parameters, config);
// Reference solid density
auto& reference_solid_density = findParameter<double>(