Potential bug in elasticity tensor?

[nilszander]

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

[koehlerson]

The discussion in kimauth/MaterialModels.jl#34 is related to this. If you call symmetric on your AD tangent, the result should be correct. Although, I can't say why this is needed

[fredrikekre]

I think that we just don't utilize in Tensors.hessian that we take the derivative of the same argument twice, so we don't know that the result should be symmetric maybe?

[nilszander]

yes! I can confirm that the test pass when computing the elasticity tensor as

∂S∂C = 4.0 * symmetric(∂²Ψ∂C²)

Many thanks!

[KristofferC]

I remember looking into this a long time ago. I don't recall the details but it had something to do with that when you take a derivative w.r.t a symmetric tensor there is some ambiguity. Googling now I find:

https://math.stackexchange.com/a/2298251

One way to interpret this is that you are thinking of 𝑇 as a function on the space of symmetric tensors. If you want to differentiate against ∂𝑇𝑘𝑙 you have to turn it into a symmetric tensor -- or possibly just a tensor. As we have seen, there are two different choices for doing so -- that is, there is some ambiguity here. And it is the details of this choice that lead to your contradiction, since depending on the choices of your identification we get one of the two answers above.

I don't fully understand the implications of all this though.

[fredrikekre]

C = transpose(F) ⋅ F

Doesn't matter for this example, but there is a special function tdot for "transpose dot product with yourself", which returns a SymmetricTensor. transpose(F) ⋅ F returns a Tensor since the ⋅ dot function do not know that the argument is the same.

[fredrikekre]

Actually, that makes this work

[nilszander]

Thanks @fredrikekre and @KristofferC.

I can confirm that using C=tdot(F) makes this work as well. To me, this suggest that the automatic differentiation works fine if we tell it that the C tensor is symmetric. However, it fails if we derive wrt a "normal" tensor, and then evaluate it on a symmetric tensor. Is this expected?

[KnutAM]

Joining the party late here, but seeing that this is still open I wanted to attempt giving a final comment that can close this issue.
Perhaps also interesting for kimauth/MaterialModels.jl#34 then?

Bonet and Wood's formulas assume all the way that C is symmetric. Doing the derivations without this assumption, I get
I = otimesl(invC', invC)

For any C, this produces that I == majortranspose(I) is true
For symmetric C, I == minortranspose(I) is true
But I[1,3,2,3] == I[1,3,3,2] is false (I[1,3,2,3] == I[3,1,3,2] is true)

Hence, I think @nilszander's last-mentioned failure is actually the result to be expected.

[nilszander]

Thanks @KnutAM for your comment. But isn't C=trans(F) . F. So to my understanding C is always symmetric. So should I see the minor and major symmetries?

[KnutAM]

When deriving the derivatives with a full tensor, no symmetry constraints are accounted for (e.g. we can perturb C_12 without touching C_21). Therefore, we get the result that I = otimesl(invC', invC) . This is equivalent to taking the derivative wrt. a C::Tensor{2} (Its numerical entries can still be symmetric)

In Bonet and Wood, however, respecting the symmetry constraints, they derived I = 0.5(otimesl(invC, invC) + otimesu(invC, invC)) This would be equivalent to taking the derivative with C::SymmetricTensor{2}

So there is, even in the theory, a difference in differentiation between being a symmetric tensor per definition and cases when a tensor just happen to have symmetric entries.

[termi-official]

A bit late to the party. The explanation on stack overflow which Kristoffer shared is very theoretic (written from a pure math stand point). Maybe a simpler explanation for people with continuum mechanics background is a paper by Itskov (see [1] below). Still very formal and not simple to read, but maybe it is more approachable to some people.

[1] Itskov, M. (2002), The Derivative with respect to a Tensor: some Theoretical Aspects and Applications. Z. angew. Math. Mech., 82: 535-544.