diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
index 938dc03730eb7d9b009804cd8461ef237d6a03af..48abe75b5852764f2440392d7d93145e55a54e1e 100644
--- a/MaterialLib/SolidModels/Ehlers-impl.h
+++ b/MaterialLib/SolidModels/Ehlers-impl.h
@@ -43,6 +43,16 @@ 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>
+ProcessLib::KelvinMatrixType<DisplacementDim> sOdotS(
+ ProcessLib::KelvinVectorType<DisplacementDim> const& v);
+
template <int DisplacementDim>
struct PhysicalStressWithInvariants final
{
@@ -216,46 +226,6 @@ void calculatePlasticResidual(
yieldFunction<DisplacementDim>(s, _mp, t, x) / G;
}
-/// Special product of \p v with itself: \f$v \odot v\f$.
-/// v is given in \p Kelvin mapping.
-template <int DisplacementDim>
-ProcessLib::KelvinMatrixType<DisplacementDim> s_odot_s(
- ProcessLib::KelvinVectorType<DisplacementDim> const& v)
-{
- ProcessLib::KelvinMatrixType<DisplacementDim> result;
-
- result(0, 0) = v(0) * v(0);
- result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
- result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
- result(0, 3) = result(3, 0) = v(0) * v(3);
- result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
- result(0, 5) = result(5, 0) = v(0) * v(5);
-
- result(1, 1) = v(1) * v(1);
- result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
- result(1, 3) = result(3, 1) = v(3) * v(1);
- result(1, 4) = result(4, 1) = v(1) * v(4);
- result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
-
- result(2, 2) = v(2) * v(2);
- result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
- result(2, 4) = result(4, 2) = v(4) * v(2);
- result(2, 5) = result(5, 2) = v(5) * v(2);
-
- result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
- result(3, 4) = result(4, 3) =
- v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
- result(3, 5) = result(5, 3) =
- v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
-
- result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
- result(4, 5) = result(5, 4) =
- v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
-
- result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
- return result;
-}
-
template <int DisplacementDim>
void calculatePlasticJacobian(
double const dt,
@@ -351,7 +321,7 @@ void calculatePlasticJacobian(
one_gt.pow_m_p * (P_dev + s.D * gm_p * M0.transpose());
// second derivative of theta
KelvinMatrix const d2theta_dsigma2 =
- theta * P_dev * s_odot_s<DisplacementDim>(sigma_D_inverse) * P_dev +
+ 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() +
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
index af91297bbd6d64e3b9dacfdc46027dd105d0e87b..3151990fadbadffe59fec5b215036e3baf029383 100644
--- a/MaterialLib/SolidModels/Ehlers.cpp
+++ b/MaterialLib/SolidModels/Ehlers.cpp
@@ -18,6 +18,67 @@ namespace Ehlers
template class SolidEhlers<2>;
template class SolidEhlers<3>;
+template <>
+ProcessLib::KelvinMatrixType<3> sOdotS<3>(
+ ProcessLib::KelvinVectorType<3> const& v)
+{
+ ProcessLib::KelvinMatrixType<3> result;
+
+ result(0, 0) = v(0) * v(0);
+ result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
+ result(0, 2) = result(2, 0) = v(5) * v(5) / 2.;
+ result(0, 3) = result(3, 0) = v(0) * v(3);
+ result(0, 4) = result(4, 0) = v(3) * v(5) / std::sqrt(2.);
+ result(0, 5) = result(5, 0) = v(0) * v(5);
+
+ result(1, 1) = v(1) * v(1);
+ result(1, 2) = result(2, 1) = v(4) * v(4) / 2.;
+ result(1, 3) = result(3, 1) = v(3) * v(1);
+ result(1, 4) = result(4, 1) = v(1) * v(4);
+ result(1, 5) = result(5, 1) = v(3) * v(4) / std::sqrt(2.);
+
+ result(2, 2) = v(2) * v(2);
+ result(2, 3) = result(3, 2) = v(5) * v(4) / std::sqrt(2.);
+ result(2, 4) = result(4, 2) = v(4) * v(2);
+ result(2, 5) = result(5, 2) = v(5) * v(2);
+
+ result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
+ result(3, 4) = result(4, 3) =
+ v(3) * v(4) / 2. + v(5) * v(1) / std::sqrt(2.);
+ result(3, 5) = result(5, 3) =
+ v(0) * v(4) / std::sqrt(2.) + v(3) * v(5) / 2.;
+
+ result(4, 4) = v(1) * v(2) + v(4) * v(4) / 2.;
+ result(4, 5) = result(5, 4) =
+ v(3) * v(2) / std::sqrt(2.) + v(5) * v(4) / 2.;
+
+ result(5, 5) = v(0) * v(2) + v(5) * v(5) / 2.;
+ return result;
+}
+
+template <>
+ProcessLib::KelvinMatrixType<2> sOdotS<2>(
+ ProcessLib::KelvinVectorType<2> const& v)
+{
+ ProcessLib::KelvinMatrixType<2> result;
+
+ result(0, 0) = v(0) * v(0);
+ result(0, 1) = result(1, 0) = v(3) * v(3) / 2.;
+ result(0, 2) = result(2, 0) = 0;
+ result(0, 3) = result(3, 0) = v(0) * v(3);
+
+ result(1, 1) = v(1) * v(1);
+ result(1, 2) = result(2, 1) = 0;
+ result(1, 3) = result(3, 1) = v(3) * v(1);
+
+ result(2, 2) = v(2) * v(2);
+ result(2, 3) = result(3, 2) = 0;
+
+ result(3, 3) = v(0) * v(1) + v(3) * v(3) / 2.;
+
+ return result;
+}
+
} // namespace Ehlers
} // namespace Solids
} // namespace MaterialLib