diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
index 938dc03730eb7d9b009804cd8461ef237d6a03af..55ab5cd5bd896c2080b24e91b8ce7a74dac147d1 100644
--- a/MaterialLib/SolidModels/Ehlers-impl.h
+++ b/MaterialLib/SolidModels/Ehlers-impl.h
@@ -43,6 +43,14 @@ namespace Solids
 {
 namespace Ehlers
 {
+
+/// Special product of \p v with itself: \f$v \odot v\f$.
+/// v is given in \p Kelvin mapping.
+/// \note Implementation only for 2 and 3 dimensions.
+template <int DisplacementDim>
+ProcessLib::KelvinMatrixType<DisplacementDim> s_odot_s(
+    ProcessLib::KelvinVectorType<DisplacementDim> const& v);
+
 template <int DisplacementDim>
 struct PhysicalStressWithInvariants final
 {
@@ -216,46 +224,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,
diff --git a/MaterialLib/SolidModels/Ehlers.cpp b/MaterialLib/SolidModels/Ehlers.cpp
index af91297bbd6d64e3b9dacfdc46027dd105d0e87b..42669ce6ccb27456afd1cfc815d0e7e8ba7673b6 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> s_odot_s<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> s_odot_s<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