MarmotVoigt.h

Go to the documentation of this file.

/* ---------------------------------------------------------------------
 *                                       _
 *  _ __ ___   __ _ _ __ _ __ ___   ___ | |_
 * | '_ ` _ \ / _` | '__| '_ ` _ \ / _ \| __|
 * | | | | | | (_| | |  | | | | | | (_) | |_
 * |_| |_| |_|\__,_|_|  |_| |_| |_|\___/ \__|
 *
 * Unit of Strength of Materials and Structural Analysis
 * University of Innsbruck,
 * 2020 - today
 *
 * festigkeitslehre@uibk.ac.at
 *
 * Matthias Neuner matthias.neuner@uibk.ac.at
 * Magdalena Schreter magdalena.schreter@uibk.ac.at
 * Alexander Dummer alexander.dummer@uibk.ac.at
 *
 * This file is part of the MAteRialMOdellingToolbox (marmot).
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of marmot.
 * ---------------------------------------------------------------------
 */
 
#pragma once
#include "Marmot/MarmotJournal.h"
#include "Marmot/MarmotMath.h"
#include "Marmot/MarmotTypedefs.h"
 
#define VOIGTFROMDIM( x ) ( ( ( x * x ) + x ) >> 1 )
 
namespace Marmot {
  namespace ContinuumMechanics::VoigtNotation {
 
    /* constexpr int VoigtSize = 6; */
 
    enum VoigtSize { OneD = 1, TwoD = 3, ThreeD = 6, Axial = 4 };
 
    constexpr VoigtSize voigtSizeFromDimension( int x )
    {
      return (VoigtSize)( ( ( x * x ) + x ) >> 1 );
    }
 
    extern const Marmot::Vector6d P;
    extern const Marmot::Vector6d PInv;
 
    extern const Marmot::Vector6d I;
    extern const Marmot::Vector6d IHyd;
    extern const Matrix6d         IDev;
 
    // Plane Stress handling
 
    Eigen::Vector3d voigtToPlaneVoigt( const Marmot::Vector6d& voigt );
 
    Marmot::Vector6d planeVoigtToVoigt( const Eigen::Vector3d& voigtPlane );
 
    template < enum VoigtSize voigtSize >
    Eigen::Matrix< double, voigtSize, 1 > reduce3DVoigt( const Marmot::Vector6d& Voigt3D )
    {
      if constexpr ( voigtSize == OneD )
        return ( Eigen::Matrix< double, 1, 1 >() << Voigt3D( 0 ) ).finished();
      else if constexpr ( voigtSize == TwoD )
        return voigtToPlaneVoigt( Voigt3D );
      else if constexpr ( voigtSize == ThreeD )
        return Voigt3D;
      else
        throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << ": invalid dimension specified" );
    }
 
    template < enum VoigtSize voigtSize >
    Marmot::Vector6d make3DVoigt( const Eigen::Matrix< double, voigtSize, 1 >& Voigt )
    {
      if constexpr ( voigtSize == OneD )
        return ( Marmot::Vector6d() << Voigt( 0 ), 0, 0, 0, 0, 0 ).finished();
      else if constexpr ( voigtSize == TwoD )
        return planeVoigtToVoigt( Voigt );
      else if constexpr ( voigtSize == ThreeD )
        return Voigt;
      else
        throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << ": invalid dimension specified" );
    }
 
    // function prototypes for  Marmot::Vector6d handling
 
    template < typename T >
    Eigen::Matrix< T, 3, 3 > voigtToStrain( const Eigen::Matrix< T, 6, 1 >& voigt )
    {
      Eigen::Matrix< T, 3, 3 > strain;
      // clang-format off
      strain << voigt[0],       voigt[3] * 0.5, voigt[4] * 0.5, 
                voigt[3] * 0.5, voigt[1],       voigt[5] * 0.5, 
                voigt[4] * 0.5, voigt[5] * 0.5, voigt[2];
      // clang-format on
      return strain;
    }
 
    template < typename T >
    Eigen::Matrix< T, 3, 3 > voigtToStress( const Eigen::Matrix< T, 6, 1 >& voigt )
    {
      Eigen::Matrix< T, 3, 3 > stress;
      // clang-format off
      stress << voigt[0], voigt[3], voigt[4], 
                voigt[3], voigt[1], voigt[5], 
                voigt[4], voigt[5], voigt[2];
      // clang-format on
      return stress;
    }
 
    Marmot::Vector6d strainToVoigt( const Eigen::Matrix3d& strainTensor );
 
    template < typename T = double >
    Eigen::Matrix< T, 6, 1 > stressToVoigt( const Eigen::Matrix< T, 3, 3 >& stressTensor )
    {
      Eigen::Matrix< T, 6, 1 > stress;
      // clang-format off
            stress << stressTensor( 0, 0 ), 
                      stressTensor( 1, 1 ), 
                      stressTensor( 2, 2 ), 
                      stressTensor( 0, 1 ),
                      stressTensor( 0, 2 ), 
                      stressTensor( 1, 2 );
      // clang-format on
      return stress;
    }
 
    Eigen::Matrix< double, 6, 6 > stiffnessToVoigt( const Eigen::Tensor< double, 4 >& C );
 
    Eigen::Tensor< double, 4 > voigtToStiffness( const Eigen::Matrix< double, 6, 6 >& voigtStiffness );
 
    template < int nDim >
    Eigen::Matrix< double, nDim, nDim > stressMatrixFromVoigt(
      const Eigen::Matrix< double, VOIGTFROMDIM( nDim ), 1 >& Voigt )
    {
      if constexpr ( nDim == 1 )
        return ( Eigen::Matrix< double, nDim, nDim >() << Voigt( 0 ) ).finished();
      else if constexpr ( nDim == 2 )
        return ( Eigen::Matrix< double, nDim, nDim >() << Voigt( 0 ), Voigt( 2 ), Voigt( 2 ), Voigt( 1 ) ).finished();
      else if constexpr ( nDim == 3 )
        return voigtToStress( Voigt );
      else
        throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << ": invalid dimension specified" );
    }
 
    template < int nDim >
    Eigen::Matrix< double, VOIGTFROMDIM( nDim ), 1 > voigtFromStrainMatrix(
      const Eigen::Matrix< double, nDim, nDim >& strain )
    {
      if constexpr ( nDim == 1 )
        return ( Eigen::Matrix< double, VOIGTFROMDIM( nDim ), 1 >() << strain( 0, 0 ) ).finished();
      else if constexpr ( nDim == 2 )
        return ( Eigen::Matrix< double, VOIGTFROMDIM( nDim ), 1 >() << strain( 0, 0 ),
                 strain( 1, 1 ),
                 2 * strain( 0, 1 ) )
          .finished();
      else if constexpr ( nDim == 3 )
        return strainToVoigt( strain );
      else
        throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << ": invalid dimension specified" );
    }
 
    namespace Invariants {
 
      Eigen::Vector3d principalStrains( const Marmot::Vector6d& strain );
 
      Eigen::Vector3d principalStresses( const Marmot::Vector6d& stress );
      // principal strains calculated from haigh westergaard strains ( sorted --> e1 > e2 > e3 )
 
      Eigen::Vector3d sortedPrincipalStrains( const Marmot::Vector6d& strain );
      // principal stressDirections calculated by solving eigenvalue problem ( !NOT sorted! )
 
      Eigen::Matrix3d principalStressesDirections( const Marmot::Vector6d& stress );
 
      double vonMisesEquivalentStress( const Marmot::Vector6d& stress );
 
      double vonMisesEquivalentStrain( const Marmot::Vector6d& strain );
 
      template < typename T >
      T normStrain( const Eigen::Matrix< T, 6, 1 >& strain )
      {
        return ContinuumMechanics::VoigtNotation::voigtToStrain( strain ).norm();
      }
 
      double normStress( const Marmot::Vector6d& stress );
      // Trace of compressive strains
 
      double StrainVolumetricNegative( const Marmot::Vector6d& strain );
 
      template < typename T >
      T I1( const Eigen::Matrix< T, 6, 1 >& stress )
      {
        return stress( 0 ) + stress( 1 ) + stress( 2 );
      }
 
      double I1Strain( const Marmot::Vector6d& strain ); //
 
      template < typename T >
      T I2( const Eigen::Matrix< T, 6, 1 >& stress )
      {
        const Eigen::Matrix< T, 6, 1 >& s = stress;
        return s( 0 ) * s( 1 ) + s( 1 ) * s( 2 ) + s( 2 ) * s( 0 ) - s( 3 ) * s( 3 ) - s( 4 ) * s( 4 ) -
               s( 5 ) * s( 5 );
      }
 
      double I2Strain( const Marmot::Vector6d& strain );
 
      template < typename T >
      T I3( const Eigen::Matrix< T, 6, 1 >& stress )
      {
        const Eigen::Matrix< T, 6, 1 >& s = stress;
        return s( 0 ) * s( 1 ) * s( 2 ) + 2. * s( 3 ) * s( 4 ) * s( 5 ) - s( 0 ) * s( 5 ) * s( 5 ) -
               s( 1 ) * s( 4 ) * s( 4 ) - s( 2 ) * s( 3 ) * s( 3 );
      }
 
      double I3Strain( const Marmot::Vector6d& strain );
 
      template < typename T >
      T J2( const Eigen::Matrix< T, 6, 1 >& stress )
      {
        const T I1_ = I1( stress );
        const T I2_ = I2( stress );
        const T res = ( 1. / 3 ) * I1_ * I1_ - I2_;
        return Marmot::Math::makeReal( res ) >= 0 ? res : T( 0.0 );
      }
 
      double J2Strain( const Marmot::Vector6d& strain );
 
      template < typename T >
      T J3( const Eigen::Matrix< T, 6, 1 >& stress )
      {
        const T I1_ = I1( stress );
        const T I2_ = I2( stress );
        const T I3_ = I3( stress );
 
        return ( 2. / 27 ) * pow( I1_, 3 ) - ( 1. / 3 ) * I1_ * I2_ + I3_;
      }
 
      double J3Strain( const Marmot::Vector6d& strain );
 
      // principal values in voigt
      std::pair< Eigen::Vector3d, Eigen::Matrix< double, 3, 6 > > principalValuesAndDerivatives(
        const Eigen::Matrix< double, 6, 1 >& S );
 
    } // namespace Invariants
 
    namespace Derivatives {
 
      // derivatives of Haigh Westergaard stresses with respect to cauchy stress in eng. notation
 
      Vector6d dStressMean_dStress();
      template < typename T >
      Eigen::Matrix< T, 6, 1 > dRho_dStress( T rho, const Eigen::Matrix< T, 6, 1 >& stress )
      {
        if ( Marmot::Math::makeReal( rho ) <= 1e-16 )
          return Eigen::Matrix< T, 6, 1 >::Zero();
 
        Eigen::Matrix< T, 6, 1 > s = IDev * stress;
 
        return T( 1. / rho ) * P.array() * s.array();
      }
 
      Marmot::Vector6d dTheta_dStress( double theta, const Marmot::Vector6d& stress );
 
      double dTheta_dJ2( const Marmot::Vector6d& stress );
 
      double dTheta_dJ3( const Marmot::Vector6d& stress );
 
      double dThetaStrain_dJ2Strain( const Marmot::Vector6d& strain );
 
      double dThetaStrain_dJ3Strain( const Marmot::Vector6d& strain );
 
      Marmot::Vector6d dJ2_dStress( const Marmot::Vector6d& stress );
 
      Marmot::Vector6d dJ3_dStress( const Marmot::Vector6d& stress );
 
      Marmot::Vector6d dJ2Strain_dStrain( const Marmot::Vector6d& strain );
 
      Marmot::Vector6d dJ3Strain_dStrain( const Marmot::Vector6d& strain );
 
      Marmot::Vector6d dThetaStrain_dStrain( const Marmot::Vector6d& strain );
 
      // derivatives of principalStess with respect to stress
 
      Marmot::Matrix36 dStressPrincipals_dStress( const Marmot::Vector6d& stress );
 
      // derivatives of plastic strains with respect to strains
 
      Eigen::Vector3d dStrainVolumetricNegative_dStrainPrincipal( const Marmot::Vector6d& strain );
 
      Matrix6d dEp_dE( const Matrix6d& CelInv, const Matrix6d& Cep );
 
      RowVector6d dDeltaEpv_dE( const Matrix6d& CelInv, const Matrix6d& Cep );
 
      Marmot::Matrix36 dSortedStrainPrincipal_dStrain( const Marmot::Vector6d& dEp );
 
      RowVector6d dDeltaEpvneg_dE( const Marmot::Vector6d& dEp, const Matrix6d& CelInv, const Matrix6d& Cep );
 
    } // namespace Derivatives
 
    namespace Transformations {
 
      Matrix6d transformationMatrixStrainVoigt( const Matrix3d& transformedCoordinateSystem );
 
      Matrix6d transformationMatrixStressVoigt( const Matrix3d& transformedCoordinateSystem );
 
      Matrix36d projectVoigtStressToPlane( const Vector3d& normalVector );
 
      Matrix36d projectVoigtStrainToPlane( const Vector3d& normalVector );
 
      Marmot::Vector6d rotateVoigtStress( const Eigen::Matrix3d& Q, const Marmot::Vector6d& stress );
 
    } // namespace Transformations
 
  }   // namespace ContinuumMechanics::VoigtNotation
} // namespace Marmot