Compute eigenvalues of rank-4 symmetric tensors

[bangerth]

Symmetric rank-4 tensors are used in both elasticity and anisotropic flow laws. A common problem I've encountered is if the rank-4 tensor is no longer positive definite (i.e., the eigenvalues of the operator that maps symmetric rank-2 tensors to symmetric rank-2 tensors are not positive).

It would be quite useful if we had a function that computes the eigenvalues of these operators. For the 2d case, the following snippet works, since a 2x2x2x2 symmetric tensor can be represented as a 3x3 matrix ("Voigt notation"):

                  SymmetricTensor<2,3> ViscoTensor;

                  ViscoTensor[0][0]=anisotropic_viscosity->stress_strain_directors[q][0][0][0][0];
                  ViscoTensor[0][1]=anisotropic_viscosity->stress_strain_directors[q][0][0][1][1];
                  ViscoTensor[0][2]=anisotropic_viscosity->stress_strain_directors[q][0][0][0][1];
                  ViscoTensor[1][0]=anisotropic_viscosity->stress_strain_directors[q][1][1][0][0];
                  ViscoTensor[1][1]=anisotropic_viscosity->stress_strain_directors[q][1][1][1][1];
                  ViscoTensor[1][2]=anisotropic_viscosity->stress_strain_directors[q][1][1][0][1];
                  ViscoTensor[2][0]=anisotropic_viscosity->stress_strain_directors[q][0][1][0][0];
                  ViscoTensor[2][1]=anisotropic_viscosity->stress_strain_directors[q][0][1][1][1];
                  ViscoTensor[2][2]=anisotropic_viscosity->stress_strain_directors[q][0][1][0][1];

                  eigenvalues(ViscoTensor);​

We could just encode such a thing in a variation of the existing eigenvalues() function.

(From a discussion with @kiralyagi.)

[jppelteret]

Compaction of tensors to FullMatrices can be done with the functions in Physics::Notation::Kelvin. One just needs to see if the scaling factors for the off-diagonal components leads to the correct eigenvalues or not.

[jppelteret]

Its likely that, for this case, you also might want to investigate the acoustic tensor. See, for example, sec. 6 of

@Article{Tanaka2014a,
  author    = {Tanaka, M. and Fujikawa, M. and Balzani, D. and Schr{\"o}der, J.},
  title     = {Robust numerical calculation of tangent moduli at finite strains based on complex-step derivative approximation and its application to localization analysis},
  journal   = {Computer Methods in Applied Mechanics and Engineering},
  year      = {2014},
  volume    = {269},
  pages     = {454--470},
  doi       = {10.1016/j.cma.2013.11.005},
}

and

@Article{Miehe2015a,
  author    = {Miehe, C. and Vallicotti, D. and Z{\"a}h, D.},
  title     = {Computational structural and material stability analysis in finite electro-elasto-statics of electro-active materials},
  journal   = {International Journal for Numerical Methods in Engineering},
  year      = {2015},
  volume    = {102},
  number    = {10},
  pages     = {1605--1637},
  doi       = {10.1002/nme.4855},
}

as well as papers by Ogden, Destrade, & co. on the same topic.

[bangerth]

I thought this was going to be easy for 2x2x2x2 symmetric tensors, given that they can be represented as a 3x3 matrix. But curiously, this 3x3 matrix may not be symmetric, so the eigenvalues are not necessarily real. We have a way to get the real-valued type for an arbitrary number type, but do we have a way to get the complex-valued type for an arbitrary type?