Matrix assembly

\bar d is likely a vector of size n2. I can’t make sense of that expression otherwise.

Then the meaning is rather clear: I’ll write it for you with more indices and summation over repeated indices is implied (I also call N_d just N so it’s clearer what the indices are):

\frac{\partial\sigma_{i,j}}{\partial \bar d_k} = -2(1- N_l \bar d_l) \sigma_{i,j} N_k

This object indeed has 3 indices and is thus a 3-tensor.

I am not sure what you mean with “implementing it in Julia” but Julia supports multi-dimensional arrays out of the box. So you could use just a 3 dimensional array to store all the values.