Potential bug in elasticity tensor?

Dear developers,
many thanks for this great tool! I just started trying out ferrite and I am impressed by the scope and the simplicity. However, I just stumbled over one question, which might be a potential bug in the computation of the rank-4 elasticity tensor.

I wrote a simple unit test to validate the correctness and used the book of Bonet and Wood as reference (see code below). While for most entries (I,J,K,L) the results match my reference solutions, the ferrite results deviate for certain indices. In particular, the ferrite results seems to be not fully symmetric as they should be according to BonetWood, Equation 6.37. As an example, consider the entry `∂S∂C[1,3,2,3]`. Due to symmetry, this entry should match `∂S∂C[1,3,3,2]`. However, in the example stated below, this is not the case. Instead, the results read
```
∂S∂C[1,3,2,3] = 0
∂S∂C[1,3,3,2] = -30.258
```
whereas the reference solution `D` reads
```
D[1,3,2,3] = -15.129
D[1,3,3,2] = -15.129
```
So it seems as if ferrite "lumps" the two entries to one.

While I don't know whether this will actually cause a problem for hyper-elasticity, it is certainly surprising. Could you help me to understand whether I am using the code correctly?

Many thanks
Nils

```
using Tensors
using Ferrite
using Test

using Tensors: Tensor

# All references refere to
# Bonet, Wood, 2016, Nonlinear Solid Mechanics for Finite Element Analysis: Statics

struct NeoHooke
μ::Float64
λ::Float64
end

function strainEnergy_NH_logJ(C, mp::NeoHooke)
μ = mp.μ
λ = mp.λ
Ic = tr(C)
J = sqrt(det(C))
return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2 # BonetWood, Eq. 6.27
end

@testset "Potential Bug" begin

μ = 123
λ = 456
material = NeoHooke(μ,λ)

# Test setup: pure shear straining ------------------------------------------------------------
g = 0.123
F = Tensor{2, 3}([1, 0, 0,
g, 1, 0,
0, 0, 1])
J = 1

# analytical solution -------------------------------------------------------------------------
C = transpose(F) ⋅ F # Equation (4.15)
invC = inv(C)
I = Tensor{4,3}((I,J,K,L) -> 0.5*(invC[I,K]*invC[J,L] + invC[I,L]*invC[J,K])) # Equation (6.36)
D = λ * invC ⊗ invC + 2*(μ-λ*log(J))*I # Equation (6.30)

# numerial solution --------------------------------------------------------------------------
strainEnergyPotential = c -> strainEnergy_NH_logJ(c, material)
∂²Ψ∂C², ∂Ψ∂C = Tensors.hessian(y -> strainEnergyPotential(y), C, :all)
S = 2.0 * ∂Ψ∂C
∂S∂C = 4.0 * ∂²Ψ∂C²

# assert ---------------------------------------------------------------------------------------
@test D[1,3,2,3] ≈ D[1,3,3,2] # symmetry of reference solution: works
@test ∂S∂C[1,3,2,3] ≈ ∂S∂C[1,3,3,2] # symmetry of numerial solution: fails
@test ∂S∂C[1,3,2,3] ≈ D[1,3,2,3] # comparison between numerial and reference solution: fails

end
```