Deviatoric and volumetric part of second order tensor for dimension not equal to 3

[TmNmr]

Hi,
I just stumbled over an issue which might lead to unexpected results for second order tensors with dimensions not equal to $3$.

To the best of my knowledge, the deviatoric part for a second order tensor $X \in \mathbb{R}^{d \times d}$ is defined by
$$\mathrm{dev}_d (X) = X - \frac{1}{d} \mathrm{tr} (X) I_d$$ with $I_d \in \mathbb{R}^{d\times d}$ denoting the identity tensor. My guess is that the denominator $3$ in the functions dev and mean could be changed to dim in order to make them applicable for general dimensions.

Tensors.jl/src/math_ops.jl

Lines 297 to 306 in 258540f

	@inline function dev(S::SecondOrderTensor)
	Tt = get_base(typeof(S))
	trace = tr(S) / 3
	Tt(
	@inline function(i, j)
	@inbounds v = i == j ? S[i,j] - trace : S[i,j]
	v
	end
	)
	end

Tensors.jl/src/math_ops.jl

Line 250 in 258540f

Statistics.mean(S::SecondOrderTensor) = tr(S) / 3

If interested I could work on a PR.

[termi-official]

Hi Timo! Good catch here. I agree that the deviator is only correct for dim=3 and that your suggested fix is the correct one.

Regarding the mean, I am not sure which definition is applied here. I would have expected that the mean of a second order tensor (in a staistical sense) should still be the sum of its entries over the number of entries. Am I missing something?

Edit: PRs are always welcome! Feel free to reach out if you get stuck somewhere.

[fredrikekre]

I believe these are both intentional. If you consider 2D just as a special case, with all 3-components to be zero, then it should still be divided by 3, not 2.

Regarding the mean, I am not sure which definition is applied here.

This definition comes from "mean stress" which is usually defined like this.

If interested I could work on a PR.

Would be great to clarify this in the documentation though, if you want to work on a PR!

[KnutAM]

I think quite some code depends on the current behavior, since e.g.
σ = 2G*dev(ϵ) + 3K*vol(ϵ) gives the correct answer for plane strain elasticity with the current implementation (implicitly assuming that the ϵ[3,3] exists and is zero). In the documentation, it is stated that we divide by 3: https://ferrite-fem.github.io/Tensors.jl/stable/man/other_operators/#Deviatoric-tensor.

I'm not saying that you are wrong @TmNmr regarding the definition, but changing this will silently lead to wrong results in code that currently gives the right result. So we would need to add some safeguards to avoid this if we decide to change the definition. I'm not sure how to do that?

[termi-official]

Thanks for elaborating Fredrik. Should we clarify the mean in the docs? I think this is really not obvious.

I think quite some code depends on the current behavior, since e.g. σ = 2G*dev(ϵ) + 3K*vol(ϵ) gives the correct answer for plane strain elasticity with the current implementation (implicitly assuming that the ϵ[3,3] exists and is zero). In the documentation, it is stated that we divide by 3: https://ferrite-fem.github.io/Tensors.jl/stable/man/other_operators/#Deviatoric-tensor.

I think the just stated motivation is the missing (crucual) piece of information in the doc string. Do you think we should add it here? Maybe we should also add clarification of plain strain and plane stress in Tensors.jl, since this seems to be a common point for errors when doing mechanics in the plane. Do you think we should open a separate issue on this for more discussion?

I'm not saying that you are wrong @TmNmr regarding the definition, but changing this will silently lead to wrong results in code that currently gives the right result. So we would need to add some safeguards to avoid this if we decide to change the definition. I'm not sure how to do that?

This is a good point. I guess we can still resolve this by

1. Adding a new function ddev or devd which implements exactly Timos definition.
2. Add a reference to the new function to the doc string of dev
3. Referencing dev in the doc string of the new function

What do you think @fredrikekre @KnutAM ?

[KnutAM]

I think the problem is that users would only have a problem with this if they don't expect it to be 3 by default, so I'm not sure if a new alternative name solves it. Perhaps just adding a !!! note in the docstring stating the choice of 3 is sufficient? If something else is intended, it is easy enough for the user to implement themselves efficiently.

Of course, the current behavior, that tr(dev(x))!=0 is not ideal. A long-term solution could be to deprecate dev/vol in favor of dev3 and vol3, and define deviatoric, volumetric for the definition by @TmNmr. After some time, dev and vol could error for non-3d, before finally re-introducing them as short-hands for deviatoric and volumetric after even more time. I think this would be the best behavior in the long run, but not sure if it is worth the hassle as this would require 2 breaking releases...

[koehlerson]

I'm not sure if I can follow. The current implementation is only valid if one assumes that 2D is a special case in the sense of plane strain elasticity, isnt it?

[KnutAM]

Not only: As @fredrikekre says, it just assumes that you have a 3d case with zeros for the values that are not present. So for calculating the deviatoric stress in 2d, it assumes plane stress. But in this case, dividing by dim=2 would be wrong for both plane stress and plane strain...

I'm trying to come up with an example where it is desirable to have the definition "dividing by dim", but I guess my brain is too much into mechanics to come up with a good example now. Does anyone else have a concrete example?

[KristofferC]

The current implementation is only valid if one assumes that 2D is a special case in the sense of plane strain elasticity, isnt it?

What do you mean with "valid" here? I mean, there are two questions:

• How do you define the deviatoric part of a tensor?
• What are you going to use the deviatoric part for in practice?

[koehlerson]

How do you define the deviatoric part of a tensor?

I would have thought that a (at least one) popular definition is that the trace is zero, which it isn't in the current implementation/definition combo. I know another colleague at our institute who recently struggled with the behavior that the trace isn't zero for a two dimensional deviatoric tensor. Back then I didn't know where it was coming from. To my feeling its a little arbitrary where we assume that 2D is a special case of 3D and add implicitly the third dimension and in other cases we don't, e.g. in the trace. I guess because there is some trade off in terms of generality/user friendliness.

What are you going to use the deviatoric part for in practice?

I don't know where @TmNmr needs it but I'd think that there are more continuum mechanical problems with vectorial unknowns that are beyond balance of linear momentum, maybe that's not true

I'm not saying that it needs to be changed, but I definitely think it's not an obvious behavior, especially since it isn't the same behavior everywhere in terms of adding the third dimension implicitly or not. It is stated here https://ferrite-fem.github.io/Tensors.jl/stable/man/other_operators/#Deviatoric-tensor but it could be more clear that this is the same behavior for all dimensions

[TmNmr]

My application was to study properties (notions of (semi)convexity) of parameter dependent functions $W\colon \mathbb{R}^{d \times d} \to \mathbb{R}$ like
$$W(F) = \frac{\mu}{k} e^{k \ ||\mathrm{dev}_d\log\sqrt{F^T F}||^2} $$
for dimensions two and three, where the deviatoric part was defined as mentioned above (projection onto the traceless tensors). So not included in a real application / finite element simulation (yet) and more a theoretical parameter study.
I was just a bit confused because my implementation in 2D didn't give the same results as my previous Matlab implementation, but 3D was perfectly fine. That's why I wasn't sure whether the current implementation was intentional (to be used for arbitrary dimensions) or the 2D case was just missing.