#include <ln_space.h>

Public Member Functions
	ln_space ()

void	pre_ln (Tensor< 2, 3 > &F)

void	post_ln (SymmetricTensor< 2, 3 > &stress_measure_T_sym, SymmetricTensor< 4, 3 > &elasto_plastic_tangent)

SymmetricTensor< 2, 3 >	post_transform (SymmetricTensor< 2, 3 > &ln_tensor)

SymmetricTensor< 2, 3 >	plastic_right_cauchy_green_AS (SymmetricTensor< 2, dim > plastic_hencky_strain)

template<>
void	post_ln (SymmetricTensor< 2, 3 > &stress_measure_T_sym, SymmetricTensor< 4, 3 > &elasto_plastic_tangent)

template<>
void	post_ln (SymmetricTensor< 2, 3 > &stress_measure_T_sym, SymmetricTensor< 4, 3 > &elasto_plastic_tangent)

Public Attributes
SymmetricTensor< 2, dim >	hencky_strain

SymmetricTensor< 2, 3 >	hencky_strain_3D

SymmetricTensor< 2, dim >	second_piola_stress_S

SymmetricTensor< 4, dim >	C

SymmetricTensor< 4, 3 >	C_3D

Private Attributes
Vector< double >	eigenvalues

std::vector< Tensor< 1, 3 > >	eigenvector

std::vector< SymmetricTensor< 2, 3 > >	eigenbasis

Vector< double >	ea

Vector< double >	da

Vector< double >	fa

SymmetricTensor< 4, 3 >	projection_tensor_P_sym

const double	comp_tolerance = 1e-8

Constructor & Destructor Documentation

template<int dim>

ln_space< dim >::ln_space ( )

 :
 eigenvalues(3),
 eigenvector(3),
 eigenbasis(3),
 ea(3),
 da(3),
 fa(3)
 {
 }

Member Function Documentation

template<int dim>

SymmetricTensor< 2, 3 > ln_space< dim >::plastic_right_cauchy_green_AS ( SymmetricTensor< 2, dim > plastic_hencky_strain )

References ln_space< dim >::ea, ln_space< dim >::eigenvector, and outer_product_sym().

 {
     AssertThrow( false, ExcMessage("plastic_right_cauchy_green_AS<< Sorry. Even though the algorithm "
                                    "to compute the plastic RCG tensor was tested, the implementation into this "
                                    "ln-space class has not been tested at all. So either you have enough faith to "
                                    "simply remove this AssertThrow in the code or you do some testing to validate the function yourself."));
 
     // Compute the eigenvalues and eigenvectors of the plastic hencky strain
      Vector<double> eigenvalues_pl(3);
      std::vector< Tensor<1,3> > eigenvector_pl(3);
      for (unsigned int i = 0; i < 3; ++i)
      {
         eigenvalues_pl[i] = eigenvectors(plastic_hencky_strain)[i].first;
         eigenvector_pl[i] = eigenvectors(plastic_hencky_strain)[i].second;
      }
 
     // Check if the found eigenvectors are perpendicular to each other
     // @todo Change the \a false to a criterion detecting whether the code is run in debug or release mode
      if ( false /*no debugging*/)
      {
         if ((fabs(eigenvalues_pl(0) - 1) > 1e-10)
                 && (fabs(eigenvalues_pl(1) - 1) > 1e-10)
                 && (fabs(eigenvalues_pl(2) - 1) > 1e-10))
             for (unsigned int i = 0; i < 3; ++i)
                 for (unsigned int j = i + 1; j < 3; ++j)
                     Assert( (fabs(eigenvector[i] * eigenvector[j]) < 1e-12),
                             ExcMessage("Eigenvectors are not perpendicular to each other") );
      }
 
     // Compute eigenbasis
      std::vector< SymmetricTensor<2,3> > eigenbasis_pl(3);
      for (unsigned int i = 0; i < 3; ++i)
         eigenbasis_pl[i] = outer_product_sym( eigenvector_pl[i] );
 
     // Compute the values of \a ea
      Vector<double> ea(3);
      for (unsigned int i = 0; i < 3; ++i)
         ea(i) = exp(2.0* eigenvalues_pl(i));            // diagonal
 
     // Finally compute the plastic right Cauchy-Green tensor
      SymmetricTensor<2,3> plastic_right_cauchy_green;
      for (unsigned int a = 0; a < 3; ++a)
         plastic_right_cauchy_green += ea(a) * eigenbasis_pl[a];
 
     return plastic_right_cauchy_green;
 }

template<int dim>

void ln_space< dim >::post_ln	(	SymmetricTensor< 2, 3 > &	stress_measure_T_sym,
		SymmetricTensor< 4, 3 > &	elasto_plastic_tangent
	)

template<>

void ln_space< 3 >::post_ln	(	SymmetricTensor< 2, 3 > &	stress_measure_T_sym,
		SymmetricTensor< 4, 3 > &	elasto_plastic_tangent
	)

References ln_space< dim >::C, ln_space< dim >::C_3D, ln_space< dim >::comp_tolerance, ln_space< dim >::da, ln_space< dim >::ea, ln_space< dim >::eigenbasis, ln_space< dim >::eigenvalues, ln_space< dim >::fa, get_tensor_operator_F_left(), get_tensor_operator_F_right(), get_tensor_operator_G(), outer_product_sym(), ln_space< dim >::projection_tensor_P_sym, ln_space< dim >::second_piola_stress_S, symmetrize(), and symmetry_check().

 {
     /*
      * 3. Set up coefficients \a theta, \a xi and \a eta
      */
 
     // Compute the coefficients based on the eigenvalues, eigenvectors and ea,da,fa
      Tensor<2,3> theta;
      Tensor<2,3> xi;
      double eta = 999999999.0;
 
     // For three different eigenvalues \f$ \lambda_a \neq \lambda_b \neq \lambda_c \f$
      if (
              ( !(std::fabs(eigenvalues(0) - eigenvalues(1)) < comp_tolerance) )
              &&
              ( !(std::fabs(eigenvalues(0) - eigenvalues(2)) < comp_tolerance) )
              &&
              ( !(std::fabs(eigenvalues(1) - eigenvalues(2)) < comp_tolerance) )
          )
      {
         eta = 0.0;
         for (unsigned int a = 0; a < 3; ++a)
             for (unsigned int b = 0; b < 3; ++b)
                 if (a != b)
                 {
                     theta[a][b] = (ea(a) - ea(b))
                                   / (eigenvalues(a) - eigenvalues(b));
                     xi[a][b] = (theta[a][b] - 0.5 * da(b))
                                / (eigenvalues(a) - eigenvalues(b));
 
                     for (unsigned int c = 0; c < 3; ++c)
                         if ((c != a) && (c != b))
                         {
                             eta +=
                                     ea(a)
                                     / (2.0
                                        * (eigenvalues(a)
                                           - eigenvalues(b))
                                        * (eigenvalues(a)
                                           - eigenvalues(c)));
                         }
                 }
      }
     //  For three equal eigenvalues \f$ \lambda_a = \lambda_b = \lambda_c \f$
      else if ( (std::fabs(eigenvalues(0) - eigenvalues(1)) < comp_tolerance)
                 &&
                (std::fabs(eigenvalues(1) - eigenvalues(2)) < comp_tolerance) )
      {
         eta = 0.0;
         for (unsigned int a = 0; a < 3; ++a)
         {
             for (unsigned int b = 0; b < 3; ++b)
                 if (a != b)
                 {
                     theta[a][b] = 0.5 * da(0);
                     xi[a][b] = (1.0 / 8.0) * fa(0);
                 }
         }
         eta = (1.0 / 8.0) * fa(0);
      }
 
     // For two equal eigenvalues a and b: \f$ \lambda_a = \lambda_b \neq \lambda_c \f$
      else if ( (std::fabs(eigenvalues(0) - eigenvalues(1)) < comp_tolerance)
                &&
                ( !(std::fabs(eigenvalues(1) - eigenvalues(2)) < comp_tolerance) ) )
      {
         eta = 0.0;
         for (unsigned int a = 0; a < 3; ++a)
             for (unsigned int b = 0; b < 3; ++b)
                 if ((a != b) && ((a == 2) || (b == 2)))
                 {
                     theta[a][b] = (ea(a) - ea(b))
                                   / (eigenvalues(a) - eigenvalues(b));
                     xi[a][b] = (theta[a][b] - 0.5 * da(b))
                                / (eigenvalues(a) - eigenvalues(b));
                 }
 
         theta[0][1] = 0.5 * da(0);
         theta[1][0] = theta[0][1];
         xi[0][1] = (1.0 / 8.0) * fa(0);
         xi[1][0] = xi[0][1];
         eta = xi[2][0];
      }
     // For two equal eigenvalues a and c: \f$ \lambda_a = \lambda_c \neq \lambda_b \f$
      else if ( (std::fabs(eigenvalues(0) - eigenvalues(2)) < comp_tolerance)
                &&
                (!(std::fabs(eigenvalues(1) - eigenvalues(2)) < comp_tolerance)) )
      {
         eta = 0.0;
         for (unsigned int a = 0; a < 3; ++a)
             for (unsigned int b = 0; b < 3; ++b)
                 if ( (a != b) && ((a == 1) || (b == 1)) )
                 {
                     theta[a][b] = (ea(a) - ea(b))
                                   / (eigenvalues(a) - eigenvalues(b));
                     xi[a][b] = (theta[a][b] - 0.5 * da(b))
                                / (eigenvalues(a) - eigenvalues(b));
                 }
 
         theta[0][2] = 0.5 * da(0);
         theta[2][0] = theta[0][2];
         xi[0][2] = (1.0 / 8.0) * fa(0);
         xi[2][0] = xi[0][2];
         eta = xi[1][0];
      }
     // For two equal eigenvalues b and c: \f$ \lambda_b = \lambda_c \neq \lambda_a \f$
      else if ( (std::fabs(eigenvalues(1) - eigenvalues(2)) < comp_tolerance)
                &&
                (!(std::fabs(eigenvalues(0) - eigenvalues(1)) < comp_tolerance)) )
      {
         eta = 0.0;
         for (unsigned int a = 0; a < 3; ++a)
             for (unsigned int b = 0; b < 3; ++b)
                 if ( (a != b) && ((a == 0) || (b == 0)) )
                 {
                     theta[a][b] = (ea(a) - ea(b))
                                   / (eigenvalues(a) - eigenvalues(b));
                     xi[a][b] = (theta[a][b] - 0.5 * da(b))
                                / (eigenvalues(a) - eigenvalues(b));
                 }
 
         theta[1][2] = 0.5 * da(1);
         theta[2][1] = theta[1][2];
         xi[1][2] = (1.0 / 8.0) * fa(1);
         xi[2][1] = xi[1][2];
         eta = xi[0][1];
      }
      else
      {
         deallog << "ln-space<< eigenvalues:0: " << eigenvalues[0] << std::endl;
         deallog << "ln-space<< eigenvalues:1: " << eigenvalues[1] << std::endl;
         deallog << "ln-space<< eigenvalues:2: " << eigenvalues[2] << std::endl;
         AssertThrow( false,
                      ExcMessage("ln-space<< Eigenvalue case not possible, check update_qph!") );
      }
 
     // Ensure that \a eta was initialised in one of the cases
      //AssertThrow( (eta < 9999999), ExcMessage("Eta in update_qph not initialised") );
 
 
     /*
      * 4. Lagrangian stresses and elasticity moduli
      */
 
      std::vector< SymmetricTensor<4,3> > Ma_x_Ma (3);
     // Compute projection tensor P
      Tensor<4,3> projection_tensor_P;
      for (unsigned int a = 0; a < 3; ++a)
      {
         Ma_x_Ma[a] = outer_product_sym(eigenbasis[a]);
         projection_tensor_P += da(a) * (Tensor<4,3> ) Ma_x_Ma[a];
         for (unsigned int b = 0; b < 3; ++b)
             if (b != a)
                 projection_tensor_P += theta[a][b] * get_tensor_operator_G(eigenbasis[a],eigenbasis[b]);
      }
 
     // Check whether the projecton tensor is symmetric and store it into a \a SymmetricTensor
      Assert( symmetry_check(projection_tensor_P), ExcMessage( "ln-space<< Projection tensor P is not symmetric") );
      projection_tensor_P_sym = symmetrize(projection_tensor_P);
 
     // Compute the double contraction of T and L
      Tensor<4,3> projection_tensor_T_doublecon_L;
      for (unsigned int a = 0; a < 3; ++a)
      {
         projection_tensor_T_doublecon_L += fa(a)
                                            * (stress_measure_T_sym * eigenbasis[a])
                                            * (Tensor<4,3> ) Ma_x_Ma[a];
         for (unsigned int b = 0; b < 3; ++b)
             if (b != a)
             {
                 projection_tensor_T_doublecon_L += 2.0 * xi[a][b]
                                                    * (
                                                         get_tensor_operator_F_right( eigenbasis[a], eigenbasis[b], eigenbasis[b], stress_measure_T_sym )
                                                       + get_tensor_operator_F_left(  eigenbasis[a], eigenbasis[b], eigenbasis[b], stress_measure_T_sym )
                                                       + get_tensor_operator_F_right( eigenbasis[b], eigenbasis[b], eigenbasis[a], stress_measure_T_sym )
                                                       + get_tensor_operator_F_left(  eigenbasis[b], eigenbasis[b], eigenbasis[a], stress_measure_T_sym )
                                                       + get_tensor_operator_F_right( eigenbasis[b], eigenbasis[a], eigenbasis[b], stress_measure_T_sym )
                                                       + get_tensor_operator_F_left(  eigenbasis[b], eigenbasis[a], eigenbasis[b], stress_measure_T_sym )
                                                      );
 
                 for (unsigned int c = 0; c < 3; ++c)
                     if ( (c != a) && (c != b) )
                     {
                         projection_tensor_T_doublecon_L += 2.0 * eta
                                                            * (
                                                                   get_tensor_operator_F_right( eigenbasis[a], eigenbasis[b], eigenbasis[c], stress_measure_T_sym )
                                                                 + get_tensor_operator_F_left(  eigenbasis[b], eigenbasis[c], eigenbasis[a], stress_measure_T_sym )
                                                               );
                     }
             }
      }
 
     // Check whether the tensor is symmetric and store it into a \a SymmetricTensor
      Assert( symmetry_check(projection_tensor_T_doublecon_L),
                   ExcMessage("ln-space<< Projection tensor T:L is not symmetric") );
      SymmetricTensor<4,3> projection_tensor_T_doublecon_L_sym = symmetrize(projection_tensor_T_doublecon_L);
 
     // Compute the retransformed values
      second_piola_stress_S = stress_measure_T_sym * projection_tensor_P_sym;
 
      C = projection_tensor_P_sym * elasto_plastic_tangent * projection_tensor_P_sym
                               + projection_tensor_T_doublecon_L_sym;
      C_3D = C;
 }

template<>

void ln_space< 2 >::post_ln	(	SymmetricTensor< 2, 3 > &	stress_measure_T_sym,
		SymmetricTensor< 4, 3 > &	elasto_plastic_tangent
	)

References ln_space< dim >::C, ln_space< dim >::C_3D, ln_space< dim >::comp_tolerance, ln_space< dim >::da, ln_space< dim >::ea, ln_space< dim >::eigenbasis, ln_space< dim >::eigenvalues, ln_space< dim >::fa, outer_product_sym(), ln_space< dim >::projection_tensor_P_sym, ln_space< dim >::second_piola_stress_S, symmetrize(), and symmetry_check().

 {
     /*
      * 3. Set up coefficients \a theta, \a xi and \a eta
      */
 
     // Compute the coefficients based on the eigenvalues, eigenvectors and ea,da,fa
      Vector<double> gamma (3);
      Vector<double> kappa (2);
      double beta = 999999999.0;
 
     // For two different eigenvalues \f$ \lambda_a \neq \lambda_b \f$
      if (std::fabs(eigenvalues(0) - eigenvalues(1)) > comp_tolerance)
      {
         beta = 2. * (ea[0]-ea[1]) / (eigenvalues(0) - eigenvalues(1));
         for (unsigned int a = 0; a < 3; ++a)
             gamma[a] = da[a] - beta;
         for (unsigned int alpha = 0; alpha < 2; ++alpha)
             kappa[alpha] = 2. * gamma[alpha] / (eigenvalues(0) - eigenvalues(1));
      }
     // For two equal eigenvalues a and b: \f$ \lambda_a = \lambda_b \f$
      else if (std::fabs(eigenvalues(0) - eigenvalues(1)) <= comp_tolerance)
      {
         beta = da[0];
         for (unsigned int a = 0; a < 3; ++a)
             gamma[a] = da[a] - da[0];
         for (unsigned int alpha = 0; alpha < 2; ++alpha)
             kappa[alpha] = (1.5 - (alpha+1)) * fa[0]; // index alpha starts at 0, hence (alpha+1)
      }
      else
      {
         deallog << "ln-space<< eigenvalues:0: " << eigenvalues[0] << std::endl;
         deallog << "ln-space<< eigenvalues:1: " << eigenvalues[1] << std::endl;
         deallog << "ln-space<< eigenvalues:2: " << eigenvalues[2] << std::endl;
         AssertThrow( false,
                      ExcMessage("ln-space<< Eigenvalue case not possible, check update_qph!") );
      }
 
     // Ensure that \a eta was initialised in one of the cases
      Assert( (beta < 9999999), ExcMessage("Beta in update_qph not initialised") );
 
 
     /*
      * 4. Lagrangian stresses and elasticity moduli
      */
 
     // zeta_a = T : M_a
      Vector<double> zeta_a (3);
      for ( unsigned int a=0; a<3; ++a )
          zeta_a(a) = stress_measure_T_sym * eigenbasis[a];
 
     // Compute projection tensor P
      Tensor<4,3> projection_tensor_P = beta * Tensor<4,3> (StandardTensors::II<3>());
 
      std::vector< SymmetricTensor<4,3> > Ma_x_Ma (3);
      // ToDo-optimize: check whether storing (M_a x M_a) saves some computation time, because we need it three times
      for (unsigned int a = 0; a < 3; ++a)
      {
          Ma_x_Ma[a] = outer_product_sym(eigenbasis[a]);
          projection_tensor_P += gamma[a] * (Tensor<4,3>) Ma_x_Ma[a];
      }
 
      projection_tensor_P_sym = symmetrize(projection_tensor_P);
 
     // Compute the double contraction of T and L
      // ToDo-optimize: merge the for-loops
      Tensor<4,3> projection_tensor_T_doublecon_L;
      for (unsigned int a = 0; a < 3; ++a)
      {
         projection_tensor_T_doublecon_L += fa(a)
                                            * zeta_a[a]
                                            * (Tensor<4,3>) Ma_x_Ma[a];
      }
      for (unsigned int alpha = 0; alpha < 2; ++alpha)
      {
         projection_tensor_T_doublecon_L += std::pow(-1, (alpha+1)+1)
                                            * kappa[alpha]
                                            * (
                                                zeta_a[alpha] * Tensor<4,3> (StandardTensors::II<3>())
                                                + outer_product_sym(stress_measure_T_sym,eigenbasis[alpha])
                                              );
      }
      for (unsigned int a = 0; a < 3; ++a)
          for (unsigned int alpha = 0; alpha < 2; ++alpha)
          {
              projection_tensor_T_doublecon_L
                   += std::pow(-1,alpha+1) * kappa[alpha]
                      * (
                            zeta_a[alpha] * Tensor<4,3> (Ma_x_Ma[a])
                            + zeta_a[a] * outer_product_sym(eigenbasis[alpha],eigenbasis[a])
                         );
          }
 
     // Check whether the tensor is symmetric and store it into a \a SymmetricTensor
      Assert( symmetry_check(projection_tensor_T_doublecon_L),
                   ExcMessage("ln-space<< Projection tensor T:L is not symmetric") );
      SymmetricTensor<4,3> projection_tensor_T_doublecon_L_sym = symmetrize(projection_tensor_T_doublecon_L);
 
     // Compute the retransformed values
      {
          SymmetricTensor<2,3> stress_S_3D = stress_measure_T_sym * projection_tensor_P_sym;
          second_piola_stress_S = extract_dim<2> ( stress_S_3D );
      }
      {
          C_3D = projection_tensor_P_sym * elasto_plastic_tangent * projection_tensor_P_sym
                 + projection_tensor_T_doublecon_L_sym;
 
          C = extract_dim<2> ( C_3D );
      }
 }

template<int dim>

SymmetricTensor< 2, 3 > ln_space< dim >::post_transform ( SymmetricTensor< 2, 3 > & ln_tensor )

References ln_space< dim >::projection_tensor_P_sym.

 {
     return ln_tensor * projection_tensor_P_sym;
 }

template<int dim>

void ln_space< dim >::pre_ln ( Tensor< 2, 3 > & F )

References ln_space< dim >::comp_tolerance, ln_space< dim >::da, ln_space< dim >::ea, ln_space< dim >::eigenbasis, ln_space< dim >::eigenvalues, ln_space< dim >::eigenvector, ln_space< dim >::fa, ln_space< dim >::hencky_strain, ln_space< dim >::hencky_strain_3D, outer_product_sym(), and symmetrize().

 {
     // Following
     // "Algorithms for computation of stresses and elasticity moduli in terms of Seth–Hill’s family of generalized strain tensors"
     // by Miehe&Lambrecht \n
     // Table I. Algorithm A
     /*
      * 1. Eigenvalues, eigenvalue bases and diagonal functions:
      */
 
     // Get the symmetric right cauchy green tensor
      SymmetricTensor<2,3> right_cauchy_green_sym = symmetrize( transpose(F) * F);//Physics::Elasticity::Kinematics::C(F);
 
     // Compute Eigenvalues, Eigenvectors and Eigenbasis
      {
         // Get Eigenvalues and Eigenvectors from the deal.ii function \a eigenvectors(*)
         for (unsigned int i = 0; i < 3; ++i) {
             eigenvalues[i] = eigenvectors(right_cauchy_green_sym)[i].first;
             eigenvector[i] = eigenvectors(right_cauchy_green_sym)[i].second;
         }
 
         // The deal.ii function \a eigenvectors return the EWe in descending order, but in Miehe et al. the C33 EW
         // shall belong to lambda_3, hence we search for the EW that equals C33
         // and move this to the end of the list of eigenvalues and eigenvectors
          // ToDo-optimize: if we find the EW at i=2, we can leave it there, hence the loop could be limited to (i<2)
         // @todo-optimize: Check the use of the std::swap function
         if ( dim==2 )
             for (unsigned int i = 0; i < 3; ++i)
                  if ( std::abs(eigenvalues[i]-right_cauchy_green_sym[2][2]) < comp_tolerance )
                  {
                      double tmp_EW = eigenvalues[2];
                      eigenvalues[2] = eigenvalues[i]; // truely lambda_3
                      eigenvalues[i] = tmp_EW;
 
                      Tensor<1,3> tmp_EV = eigenvector[2];
                      eigenvector[2] = eigenvector[i];
                      eigenvector[i] = tmp_EV;
                  }
 
         // Check if the found eigenvectors are perpendicular to each other
 //      if ((std::fabs(eigenvalues(0) - 1) > 1e-10)
 //          && (std::fabs(eigenvalues(1) - 1) > 1e-10)
 //          && (std::fabs(eigenvalues(2) - 1) > 1e-10)) {
 //          for (unsigned int i = 0; i < 3; ++i) {
 //              for (unsigned int j = i + 1; j < 3; ++j) {
 //                  Assert( (std::fabs(eigenvector[i] * eigenvector[j]) < 1e-12),
 //                               ExcMessage("ln-space<< Eigenvectors are not perpendicular to each other.") );
 //              }
 //          }
 //      }
 
         // Compute eigenbasis Ma: eigenbasis = eigenvector \otimes eigenvector
          for (unsigned int i = 0; i < 3; ++i)
          {
             eigenbasis[i] = outer_product_sym(eigenvector[i]);
             Assert( eigenvalues(i) >= 0.0,
                          ExcMessage("ln-space<< Eigenvalue is negativ. Check update_qph.") );
          }
      }
 
     // Compute diagonal function \a ea and its first and second derivate \a da and \a fa \n
      for (unsigned int i = 0; i < 3; ++i)
      {
         ea(i) = 0.5 * std::log( std::abs(eigenvalues(i)) ); // diagonal function
         da(i) = std::pow(eigenvalues(i), -1.0);             // first derivative of diagonal function ea
         fa(i) = -2.0 * std::pow(eigenvalues(i), -2.0);          // second derivative of diagonal function ea
         Assert( ea(i) == ea(i),
                      ExcMessage( "ln-space<< Ea is nan due to logarithm of negativ eigenvalue. Check update_qph.") );
         Assert( da(i) > 0.0,
                      ExcMessage( "ln-space<< First derivative da of diagonal function is "+std::to_string(da(i))+" < 0.0 ."
                                  "Check update_qph.") );
      }
 
     // Compute the Hencky strain
      for (unsigned int a = 0; a < 3; ++a)
         hencky_strain_3D += ea(a) * eigenbasis[a];
 
      hencky_strain = extract_dim<dim> (hencky_strain_3D);
 
     // Output-> SymmetricTensor<2, dim> hencky_strain, Vector<double> ea, da, fa,
     //          std::vector<Tensor<1, dim>> eigenvector, Vector<double> eigenvalues,
     //          std::vector< SymmetricTensor<2, dim> > eigenbasis
 }