pub fn reuss_average(c_a: &Tensor4, c_b: &Tensor4) -> Option<Tensor4>
Harmonic mean of two stiffness tensors (Reuss bound): C_Reuss = (C_a^{-1} + C_b^{-1})^{-1} / 2 (in Voigt/Kelvin form).