OpElasticTools Namespace Reference

Classes
struct	CommonData
	[Common data] More...
struct	FePrePostProcess
	Not used at this stage. Could be used to do some calculations, before assembly of local elements. More...
struct	OpCalculateLogStressTangent
struct	OpCalculateNeoHookeStressTangent
struct	OpEssentialRHS_U
struct	OpInternalForceRhs
struct	OpLogStrain
struct	OpPostProcElastic
struct	OpRotate
struct	OpSaveReactionForces
struct	OpSchurBegin
struct	OpSchurBeginBoundary
struct	OpSchurEnd
struct	OpSchurEndBoundary
struct	OpSchurPreconditionMassContact
struct	OpSchurPreconditionMassTau
struct	OpSetMaterialBlock
struct	OpSetMaterialBlockBoundary
struct	OpSetUnSetRowData
struct	OpStrain
struct	OpStress

Typedefs
using	OpMassPrecHTensorBound = FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 3, 9 >
using	OpMassPrecHTensorVol = FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 3, 9 >
using	OpMassPrecScalar = FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 1 >
typedef boost::function< Tensor1< double, 3 >(const double, const double, const double)>	VectorFun
	[Operators definitions] More...
typedef struct OpSchurPreconditionMassContact< OpMassPrecHTensorBound >	OpSchurPreconditionMassContactBound
typedef struct OpSchurPreconditionMassContact< OpMassPrecHTensorVol >	OpSchurPreconditionMassContactVol

Functions
template<typename T >
auto	to_non_symm (T &symm)
template<typename T >
auto	to_symm (T &non_symm)
bool	is_eq (const double &a, const double &b)
bool	is_diff (const double &a, const double &b)
template<typename T >
Ddg< double, 3, 3 >	tensor4_to_ddg (Tensor4< T, 3, 3, 3, 3 > &t)
template<typename T >
auto	get_tensor4_from_ddg (Ddg< T, 3, 3 > &ddg)
template<typename T >
bool	is_ddg (Tensor4< T, 3, 3, 3, 3 > tensor)
template<typename T >
auto	exp_map (const Tensor2_symmetric< T, 3 > &X)
template<typename T >
auto	get_rotation_R (Tensor1< T, 3 > t1_omega, double tt)
template<typename T >
auto	get_log_map (const Tensor2_symmetric< T, 3 > &X)
template<typename T >
auto	get_diff_log_map (Tensor2_symmetric< T, 3 > &X)
template<typename T >
MoFEMErrorCode	get_eigen_val_and_proj_lapack (const Tensor2_symmetric< T, 3 > &X, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec)
template<typename T1 , typename T2 , typename T3 >
auto	get_ln_X_lapack (const Tensor2_symmetric< T1, 3 > &X, Tensor1< T2, 3 > &lambda, Tensor2< T3, 3, 3 > &eig_vec)
template<typename T >
auto	get_ln_X (const Tensor2_symmetric< T, 3 > &X)
template<typename T , typename TP , typename T3 >
auto	calculate_nominal_moduli_TL (Tensor2_symmetric< T, 3 > &X, Tensor2< TP, 3, 3 > &F, Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Tensor4< T3, 3, 3, 3, 3 > &TLs3, bool compute_tangent)
template<typename T , typename TP , typename T3 >
auto	test_projection_lukasz (Tensor2_symmetric< T, 3 > &X, Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Ddg< T3, 3, 3 > &TLs3)
template<bool TANGENT, typename T , typename TP , typename T3 >
auto	calculate_lagrangian_moduli_TL (Tensor2_symmetric< TP, 3 > &stress, Tensor1< T, 3 > &lambda, Tensor2< T, 3, 3 > &eig_vec, Ddg< T3, 3, 3 > &TLs3)
template<typename T >
Tensor2< double, 3, 3 >	getFiniteDiff (const int &kk, const int &ll, CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)
template<typename T >
Tensor4< double, 3, 3, 3, 3 >	getDTensorFunction (CommonData &common_data, Tensor2< T, 3, 3 > &F, Tensor2< double, 3, 3 > &t_stress)

Typedef Documentation

OpMassPrecHTensorBound

using OpElasticTools::OpMassPrecHTensorBound = typedef FormsIntegrators<BoundaryEleOp>::Assembly< USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>

Definition at line 21 of file ElasticOperators.hpp.

OpMassPrecHTensorVol

using OpElasticTools::OpMassPrecHTensorVol = typedef FormsIntegrators<DomainEleOp>::Assembly< USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>

Definition at line 24 of file ElasticOperators.hpp.

OpMassPrecScalar

using OpElasticTools::OpMassPrecScalar = typedef FormsIntegrators<DomainEleOp>::Assembly< USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<1, 1>

Definition at line 27 of file ElasticOperators.hpp.

OpSchurPreconditionMassContactBound

typedef struct OpSchurPreconditionMassContact< OpMassPrecHTensorBound > OpElasticTools::OpSchurPreconditionMassContactBound

Definition at line 1194 of file ElasticOperators.hpp.

OpSchurPreconditionMassContactVol

typedef struct OpSchurPreconditionMassContact< OpMassPrecHTensorVol > OpElasticTools::OpSchurPreconditionMassContactVol

Definition at line 1195 of file ElasticOperators.hpp.

VectorFun

typedef boost::function<Tensor1<double, 3>(const double, const double, const double)> OpElasticTools::VectorFun

[Operators definitions]

Definition at line 836 of file ElasticOperators.hpp.

Function Documentation

calculate_lagrangian_moduli_TL()

template<bool TANGENT, typename T , typename TP , typename T3 >

auto OpElasticTools::calculate_lagrangian_moduli_TL ( Tensor2_symmetric< TP, 3 > &  stress, Tensor1< T, 3 > &  lambda, Tensor2< T, 3, 3 > &  eig_vec, Ddg< T3, 3, 3 > &  TLs3  )

inline

Definition at line 610 of file ElasticOperators.hpp.

                                                                {
 
  Ddg<double, 3, 3> P3;
  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
 
  VectorDouble e(3);
  VectorDouble f(3);
  VectorDouble d(3);
  MatrixDouble Vi(3, 3);
  MatrixDouble Ksi(3, 3);
  MatrixDouble Zeta(3, 3);
  Vi.clear();
  Ksi.clear();
  Zeta.clear();
  double eta = 0;
 
  Tensor1<double, 3> _N[3];
  for (int dd = 0; dd != 3; ++dd) {
    for (int cc = 0; cc != 3; ++cc)
      _N[dd](cc) = eig_vec(dd, cc);
  }
 
  for (int ii = 0; ii != 3; ++ii) {
    e(ii) = 0.5 * log(lambda(ii));
    d(ii) = 1. / lambda(ii);
    f(ii) = -2. / (lambda(ii) * lambda(ii));
  }
 
  const double &lam0 = lambda(0);
  const double &lam1 = lambda(1);
  const double &lam2 = lambda(2);
 
  auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
    Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
    Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
 
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
          if (jj == k_ || ii == k_) {
            Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
            Ksi(ii, jj) =
                (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          }
        }
    eta = Ksi(k_, i_);
  };
 
  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
    // all different
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj) {
        if (ii != jj) {
 
          Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
          Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          for (int kk = 0; kk != 3; ++kk) {
 
            if (kk != ii && kk != jj)
              eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
                              (lambda(ii) - lambda(kk)));
          }
        }
      }
 
  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
    // all the same
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
 
          Vi(ii, jj) = 0.5 * d(ii);
          Ksi(ii, jj) = 0.125 * f(ii);
        }
    eta = 0.125 * f(0);
  } else if (is_eq(lam0, lam1)) {
 
    compute_ksi_eta_2eq(0, 1, 2);
 
  } else if (is_eq(lam0, lam2)) {
 
    compute_ksi_eta_2eq(0, 2, 1);
 
  } else if (is_eq(lam1, lam2)) {
 
    compute_ksi_eta_2eq(1, 2, 0);
  }
 
  P3(i, j, k, l) = 0.;
  Tensor2_symmetric<double, 3> Mii;
  Tensor2_symmetric<double, 3> Mij;
  Tensor2_symmetric<double, 3> Mik;
  Tensor2_symmetric<double, 3> Mjk;
  Tensor2_symmetric<double, 3> Mjj;
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
    P3(i, j, k, l) += ((0.5 * d(ii) * _N[ii](i)) ^ _N[ii](j)) * Mii(k, l);
    for (int jj = 0; jj != 3; ++jj)
      if (ii != jj) {
        Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
        P3(i, j, k, l) += ((Vi(ii, jj) * _N[ii](i)) ^ _N[jj](j)) * Mij(k, l);
      }
  }
 
  if (!TANGENT)
    return P3;
 
  TLs3(i, j, k, l) = 0.;
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      Zeta(ii, jj) = stress(i, j) * (_N[ii](i) ^ _N[jj](j));
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
    TLs3(i, j, k, l) += (0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j)) * Mii(k, l);
  }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii) {
        Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
        Mjj(i, j) = (_N[jj](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[jj](j));
        const double k2 = 2. * Ksi(ii, jj);
        const double zp = k2 * Zeta(ii, jj);
        const double zp2 = k2 * Zeta(jj, jj);
        TLs3(i, j, k, l) += (zp * Mij(i, j)) * Mjj(k, l) +
                            (zp * Mjj(i, j)) * Mij(k, l) +
                            (zp2 * Mij(i, j)) * Mij(k, l);
        for (int kk = 0; kk != 3; ++kk) {
          if (kk != ii && kk != jj) {
            Mik(i, j) = (_N[ii](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[ii](j));
            Mjk(i, j) = (_N[jj](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[jj](j));
            TLs3(i, j, k, l) +=
                (2. * eta * Zeta(ii, jj) * Mik(i, j)) * Mjk(k, l);
          }
        }
      }
 
  return P3;
};

calculate_nominal_moduli_TL()

template<typename T , typename TP , typename T3 >

auto OpElasticTools::calculate_nominal_moduli_TL ( Tensor2_symmetric< T, 3 > &  X, Tensor2< TP, 3, 3 > &  F, Tensor2_symmetric< TP, 3 > &  stress, Tensor1< T, 3 > &  lambda, Tensor2< T, 3, 3 > &  eig_vec, Tensor4< T3, 3, 3, 3, 3 > &  TLs3, bool  compute_tangent  )

inline

Definition at line 382 of file ElasticOperators.hpp.

                                                  {
 
  // Tensor4<double, 3, 3, 3, 3> TLs3;
  Tensor4<double, 3, 3, 3, 3> P3;
  Tensor2<double, 3, 3> invF;
 
  double det_f;
  CHKERR determinantTensor3by3(F, det_f);
  CHKERR invertTensor3by3(F, det_f, invF);
 
  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
 
  // VectorDouble lambda(3);
  VectorDouble e(3);
  VectorDouble f(3);
  VectorDouble d(3);
  MatrixDouble Vi(3, 3);
  MatrixDouble Ksi(3, 3);
  MatrixDouble Zeta(3, 3);
  Vi.clear();
  Ksi.clear();
  Zeta.clear();
  double eta = 0;
 
  Tensor1<double, 3> _N[3];
  Tensor1<double, 3> _n[3];
 
  for (int dd : {0, 1, 2}) {
    for (int cc : {0, 1, 2})
      _N[dd](cc) = eig_vec(dd, cc);
  }
 
  for (int dd : {0, 1, 2})
    _n[dd](i) = F(i, j) * _N[dd](j);
 
  for (int ii = 0; ii != 3; ++ii) {
    e(ii) = 0.5 * log(lambda(ii));
    d(ii) = 1. / lambda(ii);
    f(ii) = -2. / (lambda(ii) * lambda(ii));
  }
 
  const double &lam0 = lambda(0);
  const double &lam1 = lambda(1);
  const double &lam2 = lambda(2);
 
  auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
    Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
    Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
 
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
          if (jj == k_ || ii == k_) {
            Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
            Ksi(ii, jj) =
                (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          }
        }
    eta = Ksi(k_, i_);
  };
 
  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
    // all different
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj) {
        if (ii != jj) {
 
          Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
          Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
          for (int kk = 0; kk != 3; ++kk) {
 
            if (kk != ii && kk != jj)
              eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
                              (lambda(ii) - lambda(kk)));
          }
        }
      }
 
  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
    // all the same
    for (int ii = 0; ii != 3; ++ii)
      for (int jj = 0; jj != 3; ++jj)
        if (ii != jj) {
 
          Vi(ii, jj) = 0.5 * d(ii);
          Ksi(ii, jj) = 0.125 * f(ii);
        }
    eta = 0.125 * f(0);
  } else if (is_eq(lam0, lam1)) {
 
    compute_ksi_eta_2eq(0, 1, 2);
 
  } else if (is_eq(lam0, lam2)) {
 
    compute_ksi_eta_2eq(0, 2, 1);
 
  } else if (is_eq(lam1, lam2)) {
 
    compute_ksi_eta_2eq(1, 2, 0);
  }
 
  P3(i, j, k, l) = 0.;
 
  Tensor2<double, 3, 3> Mii;
  Tensor2<double, 3, 3> Mij;
  Tensor2<double, 3, 3> Mik;
  Tensor2<double, 3, 3> Mjk;
  Tensor2<double, 3, 3> Mjj;
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
    P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
    for (int jj = 0; jj != 3; ++jj)
      if (ii != jj) {
        Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
        P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
      }
  }
 
  if (!compute_tangent)
    return P3;
 
  TLs3(i, j, k, l) = 0.;
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
 
  for (int ii = 0; ii != 3; ++ii) {
    Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
    TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
  }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii)
        for (int kk = 0; kk != 3; ++kk)
          if (kk != ii && kk != jj) {
            Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
            Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
            TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
          }
 
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      if (jj != ii) {
        Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
        Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
        const double k2 = 2. * Ksi(ii, jj);
        const double zp = k2 * Zeta(ii, jj);
        const double zp2 = k2 * Zeta(jj, jj);
        TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
                            zp * Mjj(i, j) * Mij(k, l) +
                            zp2 * Mij(i, j) * Mij(k, l);
      }
 
  Tensor2<double, 3, 3> Piola;
  Tensor2<double, 3, 3> invFPiola;
  Tensor4<double, 3, 3, 3, 3> Ttmp3;
  auto nstress = to_non_symm(stress);
  Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
  invFPiola(i, j) = invF(i, k) * Piola(k, j);
  Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
                      0.5 * invFPiola(i, l) * kronecker_delta(j, k);
  TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
 
  return P3;
};

exp_map()

template<typename T >

auto OpElasticTools::exp_map ( const Tensor2_symmetric< T, 3 > &  X )

inline

Definition at line 200 of file ElasticOperators.hpp.

                                                                            {
  Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
  Tensor2<double, 3, 3> X_nt;
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  int fac_n = 1;
 
  auto Xcopy = to_non_symm(X);
  auto X_n = to_non_symm(X);
  exp_tensor(i, j) += Xcopy(i, j);
 
  for (int n = 2; n != max_it; ++n) {
    fac_n *= n;
    X_nt(i, j) = X_n(i, j);
    X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
    const double inv_nf = 1. / fac_n;
    exp_tensor(i, j) += inv_nf * X_n(i, j);
    const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
    if (e < eps)
      return to_symm(exp_tensor);
  }
 
  return to_symm(exp_tensor);
};

get_diff_log_map()

template<typename T >

auto OpElasticTools::get_diff_log_map ( Tensor2_symmetric< T, 3 > &  X )

inline

Definition at line 285 of file ElasticOperators.hpp.

                                                                               {
  Tensor2<double, 3, 3> X_nt_tmp;
  Tensor4<double, 3, 3, 3, 3> diff_log_map;
  Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  auto get_Xn = [&](auto &log_tensor) {
    vector<Tensor2<double, 3, 3>> Xn;
    Xn.emplace_back();
    auto &cu_Xn = Xn.back();
    cu_Xn(i, j) = kronecker_delta(i, j);
    for (int n = 2; n != max_it; ++n) {
      const double cof = 1. / n;
      auto X_nt_tmp = Xn.back();
      Xn.emplace_back();
      auto &c_Xn = Xn.back();
      c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
      const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
      if (e < eps)
        return Xn;
    }
    return Xn;
  };
 
  auto log_tensor = to_non_symm(X);
  log_tensor(i, j) -= kronecker_delta(i, j);
  diff_log_map(i, j, k, l) = 0.;
 
  auto Xn_vec = get_Xn(log_tensor);
 
  for (int n = 1; n != Xn_vec.size() + 1; ++n) {
    const double cof = pow(-1, n + 1) / n;
    for (int m = 1; m != n + 1; ++m) {
      auto &Xn_1 = Xn_vec[m - 1];
      auto &Xn_2 = Xn_vec[n - m];
      diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
    }
  }
 
  auto isddg = tensor4_to_ddg((*cache).Is);
  diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
  return tensor4_to_ddg(diff_log_map_tmp);
};

get_eigen_val_and_proj_lapack()

template<typename T >

MoFEMErrorCode OpElasticTools::get_eigen_val_and_proj_lapack ( const Tensor2_symmetric< T, 3 > &  X, Tensor1< T, 3 > &  lambda, Tensor2< T, 3, 3 > &  eig_vec  )

inline

Definition at line 333 of file ElasticOperators.hpp.

                                                         {
 
  MoFEMFunctionBeginHot;
  Tensor1<double, 3> eig;
  auto eigen_vec_lap = to_non_symm(X);
  constexpr int n = 3;
  constexpr int lda = 3;
  constexpr int lwork = 15;
  std::array<double, 15> work;
  if (lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0),
                   work.data(), lwork) > 0)
    SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
            "The algorithm failed to compute eigenvalues.");
 
  lambda(i) = eig(i);
  eig_vec(i, j) = eigen_vec_lap(i, j);
 
  // perturb
  // const double h = 1e-12;
  // if (is_eq(lambda[0], lambda[1]))
  //   lambda[0] *= (1 + h);
  // if (is_eq(lambda[1], lambda[2]))
  //   lambda[2] *= (1 - h);
 
  MoFEMFunctionReturnHot(0);
}

get_ln_X()

template<typename T >

auto OpElasticTools::get_ln_X ( const Tensor2_symmetric< T, 3 > &  X )

inline

Definition at line 376 of file ElasticOperators.hpp.

                                                                             {
  return get_log_map(X);
};

get_ln_X_lapack()

template<typename T1 , typename T2 , typename T3 >

auto OpElasticTools::get_ln_X_lapack ( const Tensor2_symmetric< T1, 3 > &  X, Tensor1< T2, 3 > &  lambda, Tensor2< T3, 3, 3 > &  eig_vec  )

inline

Definition at line 363 of file ElasticOperators.hpp.

                                                        {
  Tensor2_symmetric<double, 3> lnX;
  Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
  for (int ii = 0; ii != 3; ++ii)
    diag(ii, ii) = log(lambda(ii));
 
  lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
 
  return lnX;
};

get_log_map()

template<typename T >

auto OpElasticTools::get_log_map ( const Tensor2_symmetric< T, 3 > &  X )

inline

Definition at line 259 of file ElasticOperators.hpp.

                                                          {
 
  constexpr int max_it = 1000;
  constexpr double eps = 1e-12;
 
  auto log_tensor = to_non_symm(X);
  log_tensor(i, j) -= kronecker_delta(i, j);
 
  auto Xcopy = to_non_symm(log_tensor);
  auto X_n = to_non_symm(log_tensor);
  Tensor2<double, 3, 3> X_nt_tmp;
 
  for (int n = 2; n != max_it; ++n) {
 
    X_nt_tmp(i, j) = X_n(i, j);
    X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
    const double inv_n = 1. / n;
    log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
    const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
    if (e < eps)
      break;
  }
 
  return to_symm(log_tensor);
};

get_rotation_R()

template<typename T >

auto OpElasticTools::get_rotation_R ( Tensor1< T, 3 >  t1_omega, double  tt  )

inline

Definition at line 228 of file ElasticOperators.hpp.

                                                              {
 
  auto get_rotation = [&](Tensor1<double, 3> &t_omega) {
    FTensor::Tensor2<double, 3, 3> t_R;
 
    constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
    t_R(i, j) = t_kd(i, j);
 
    const double angle = sqrt(t_omega(i) * t_omega(i));
    if (std::abs(angle) < 1e-18)
      return t_R;
 
    FTensor::Tensor2<double, 3, 3> t_Omega;
    t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
    const double a = sin(angle) / angle;
    const double ss_2 = sin(angle / 2.);
    const double b = 2. * ss_2 * ss_2 / (angle * angle);
    t_R(i, j) += a * t_Omega(i, j);
    t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
 
    return t_R;
  };
 
  Tensor1<double, 3> my_c_omega;
  my_c_omega(i) = t1_omega(i) * tt;
  auto t_R = get_rotation(my_c_omega);
 
  return get_rotation(my_c_omega);
};

get_tensor4_from_ddg()

template<typename T >

auto OpElasticTools::get_tensor4_from_ddg ( Ddg< T, 3, 3 > &  ddg )

inline

Definition at line 175 of file ElasticOperators.hpp.

                                                                          {
  Tensor4<double, 3, 3, 3, 3> tens4;
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      for (int kk = 0; kk != 3; ++kk)
        for (int ll = 0; ll != 3; ++ll)
          tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
  return tens4;
};

getDTensorFunction()

template<typename T >

Tensor4< double, 3, 3, 3, 3 > OpElasticTools::getDTensorFunction ( CommonData &  common_data, Tensor2< T, 3, 3 > &  F, Tensor2< double, 3, 3 > &  t_stress  )

inline

Definition at line 817 of file ElasticOperators.hpp.

                                                    {
  Tensor4<double, 3, 3, 3, 3> my_D;
  for (int k = 0; k != 3; ++k) {
    for (int l = 0; l != 3; ++l) {
      auto d_stress = getFiniteDiff(k, l, common_data, F, t_stress);
      for (int i = 0; i != 3; ++i) {
        for (int j = 0; j != 3; ++j) {
          my_D(i, j, k, l) = d_stress(i, j);
        }
      }
    }
  }
  return my_D;
};

getFiniteDiff()

template<typename T >

Tensor2< double, 3, 3 > OpElasticTools::getFiniteDiff ( const int &  kk, const int &  ll, CommonData &  common_data, Tensor2< T, 3, 3 > &  F, Tensor2< double, 3, 3 > &  t_stress  )

inline

Definition at line 757 of file ElasticOperators.hpp.

                                                                    {
  MoFEMFunctionBeginHot;
  Tensor2<double, 3, 3> out;
 
  Tensor2<double, 3, 3> strainPl;
  Tensor2<double, 3, 3> strainMin;
  strainPl(i, j) = F(i, j);
  strainMin(i, j) = F(i, j);
  const double h = 1e-8;
  strainPl(kk, ll) += h;
  strainMin(kk, ll) -= h;
  Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
 
  auto calculate_stress = [&](auto &F) {
    Tensor2<double, 3, 3> Piola;
    Tensor2_symmetric<double, 3> stress_sym;
    Tensor2_symmetric<double, 3> C;
    Tensor2_symmetric<double, 3> strain;
 
    C(i, j) = F(m, i) ^ F(m, j);
 
    Tensor1<double, 3> lambda;
    Tensor2<double, 3, 3> eig_vec;
    CHKERR get_eigen_val_and_proj_lapack(C, lambda, eig_vec);
 
    auto lnC = get_ln_X_lapack(C, lambda, eig_vec);
    // auto lnC = get_log_map(C);
    strain(i, j) = 0.5 * lnC(i, j);
    auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*common_data.mtD);
    stress_sym(i, j) = t_D(i, j, k, l) * strain(k, l);
    auto stress_tmp = to_non_symm(stress_sym);
 
    // Tensor4<double, 3, 3, 3, 3> D2;
    // constexpr auto t_kd = Kronecker_Delta<int>();
    // Tensor4<double, 3, 3, 3, 3> dC_dF;
    // dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
    // auto t_L = get_diff_log_map(C);
    // t_L(i, j, k, l) *= 0.5;
    // D2(i, j, k, l) = t_L(i, j, m, n) * dC_dF(m, n, k, l);
 
    auto D2 = calculate_nominal_moduli_TL(C, F, stress_sym, lambda, eig_vec,
                                          dummy, false);
 
    Piola(k, l) = stress_tmp(i, j) * D2(i, j, k, l);
 
    return Piola;
  };
 
  auto StressPlusH = calculate_stress(strainPl);
  auto StressMinusH = calculate_stress(strainMin);
  // out(i, j) = (StressPlusH(i, j) - t_stress(i, j)) / (h);
 
  out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
 
  return out;
}

is_ddg()

template<typename T >

bool OpElasticTools::is_ddg ( Tensor4< T, 3, 3, 3, 3 >  tensor )

inline

Definition at line 185 of file ElasticOperators.hpp.

                                                                        {
  auto test_ddg = tensor4_to_ddg(tensor);
  for (int ii = 0; ii != 3; ++ii)
    for (int jj = 0; jj != 3; ++jj)
      for (int kk = 0; kk != 3; ++kk)
        for (int ll : {0, 1, 2}) {
          const double diff =
              abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
          if (diff > 1e-6)
            return false;
        }
 
  return true;
};

is_diff()

bool OpElasticTools::is_diff ( const double &  a, const double &  b  )

inline

Definition at line 122 of file ElasticOperators.hpp.

                                                      {
  constexpr double eps = 1e-12;
  return abs(a - b) > eps;
};

is_eq()

bool OpElasticTools::is_eq ( const double &  a, const double &  b  )

inline

Definition at line 118 of file ElasticOperators.hpp.

                                                    {
  constexpr double eps = 1e-12;
  return abs(a - b) < eps;
};

tensor4_to_ddg()

template<typename T >

Ddg< double, 3, 3 > OpElasticTools::tensor4_to_ddg ( Tensor4< T, 3, 3, 3, 3 > &  t )

inline

Definition at line 128 of file ElasticOperators.hpp.

                                                                   {
  Ddg<double, 3, 3> loc_tens;
 
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
 
  loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
  loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
  loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
  loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
  loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
  loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
  loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
  loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
  loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
  loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
  loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
  loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
  loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
  loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
  loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
  loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
  loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
  loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
  loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
  loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
  loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
  loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
  loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
  loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
  loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
  loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
  loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
  loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
  loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
  loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
  loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
  loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
  loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
  loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
  loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
  loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
 
  return loc_tens;
};

test_projection_lukasz()

template<typename T , typename TP , typename T3 >

auto OpElasticTools::test_projection_lukasz ( Tensor2_symmetric< T, 3 > &  X, Tensor2_symmetric< TP, 3 > &  stress, Tensor1< T, 3 > &  lambda, Tensor2< T, 3, 3 > &  eig_vec, Ddg< T3, 3, 3 > &  TLs3  )

inline

Definition at line 556 of file ElasticOperators.hpp.

                                                                       {
 
  Tensor1<double, 3> eig{lambda(0), lambda(1), lambda(2)};
  Tensor2<double, 3, 3> eigen_vec;
  eigen_vec(i, j) = eig_vec(i, j);
 
  auto get_uniq_nb = [&](auto eig) {
    return distance(&eig(0), unique(&eig(0), &eig(0) + 3, is_eq));
  };
 
  auto sort_eigen_vals2 = [&](auto &eig, auto &eigen_vec) {
    if (is_eq(eig(0), eig(1))) {
      Tensor2<double, 3, 3> eigen_vec_c{
          eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
          eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
          eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
      Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
      eig(i) = eig_c(i);
      eigen_vec(i, j) = eigen_vec_c(i, j);
    } else {
      Tensor2<double, 3, 3> eigen_vec_c{
          eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
          eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
          eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
      Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
      eig(i) = eig_c(i);
      eigen_vec(i, j) = eigen_vec_c(i, j);
    }
  };
 
  Tensor2_symmetric<double, 3> Ts;
  Ts(i, j) = stress(i, j);
  auto f = [](double v) { return 0.5 * log(v); };
  auto d_f = [](double v) { return 0.5 / v; };
  auto dd_f = [](double v) { return -0.5 / (v * v); };
 
  auto nb_uniq = get_uniq_nb(eig);
  if (nb_uniq == 2)
    sort_eigen_vals2(eig, eigen_vec);
  auto t_d =
 
      EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
 
  auto t_dd =
 
      EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, Ts, nb_uniq);
 
  TLs3(i, j, k, l) = t_dd(i, j, k, l);
  return t_d;
};

to_non_symm()

template<typename T >

auto OpElasticTools::to_non_symm ( T &  symm )

inline

Definition at line 90 of file ElasticOperators.hpp.

                                                       {
  Tensor2<double, 3, 3> non_symm;
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
  non_symm(N0, N0) = symm(N0, N0);
  non_symm(N1, N1) = symm(N1, N1);
  non_symm(N2, N2) = symm(N2, N2);
  non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
  non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
  non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
  return non_symm;
};

to_symm()

template<typename T >

auto OpElasticTools::to_symm ( T &  non_symm )

inline

Definition at line 104 of file ElasticOperators.hpp.

                                                       {
  Tensor2_symmetric<double, 3> symm;
  Number<0> N0;
  Number<1> N1;
  Number<2> N2;
  symm(N0, N0) = non_symm(N0, N0);
  symm(N1, N1) = non_symm(N1, N1);
  symm(N2, N2) = non_symm(N2, N2);
  symm(N0, N1) = non_symm(N0, N1);
  symm(N0, N2) = non_symm(N0, N2);
  symm(N1, N2) = non_symm(N1, N2);
  return symm;
};