materials: build the deviatoric projector once, not per call (#77)

include/numsim-materials/materials/drucker_prager_yield_function.h

@@ -157,7 +157,7 @@ struct drucker_prager_yield_function {
   /// Takes sig_dev directly (not N) to avoid cancellation error from
   /// reconstructing s = (N - β/3·I)·2q when q is small.
   tensor4 flow_normal_stress_derivative(const tensor2& sig_dev, T sqrt_j2) const {
-    const tensor4 IIdev{plasticity_detail::make_IIdev<T, Dim>()};
+    const auto& IIdev = plasticity_detail::cached_IIdev<T, Dim>();
     const auto j2 = sqrt_j2 * sqrt_j2;
     return (IIdev - tmech::otimes(sig_dev, sig_dev) / (T{2} * j2)) / (T{2} * sqrt_j2);
   }

include/numsim-materials/materials/plasticity_utils.h

@@ -17,6 +17,13 @@ tmech::tensor<T, Dim, 4> make_IIdev() {
   return tmech::tensor<T, Dim, 4>{IIsym - IIvol};
 }
 
+/// The same projector, built once per (T, Dim) and shared.
+template<typename T, std::size_t Dim>
+const tmech::tensor<T, Dim, 4>& cached_IIdev() {
+  static const tmech::tensor<T, Dim, 4> value{make_IIdev<T, Dim>()};
+  return value;
+}
+
 /// Isotropic elastic tangent C = K*I⊗I + 2G*IIdev.
 ///
 /// The bulk term is K*I⊗I, NOT 3K*IIvol. Those are equal only at Dim = 3,

include/numsim-materials/materials/yield_functions.h

@@ -64,7 +64,7 @@ struct j2_yield_function {
   /// Takes sig_dev (not N) to match the DP interface and avoid
   /// needing to reconstruct s from N.
   tensor4 flow_normal_stress_derivative(const tensor2& sig_dev, T sig_eq) const {
-    const tensor4 IIdev{plasticity_detail::make_IIdev<T, Dim>()};
+    const auto& IIdev = plasticity_detail::cached_IIdev<T, Dim>();
     const tensor2 N{T{1.5} * sig_dev / sig_eq};
     return (T{1.5} * IIdev - tmech::otimes(N, N)) / sig_eq;
   }