v0.10.0
Classes | Namespaces | Macros | Functions
ElasticOperators.hpp File Reference

Go to the source code of this file.

Classes

struct  OpElasticTools::CommonData
 [Common data] More...
 
struct  OpElasticTools::OpStrain
 
struct  OpElasticTools::OpLogStrain
 
struct  OpElasticTools::OpStress
 
struct  OpElasticTools::OpCalculateLogStressTangent
 
struct  OpElasticTools::OpRotatingFrameRhs
 
struct  OpElasticTools::OpRotatingFrameBoundaryRhs
 
struct  OpElasticTools::OpRotatingFrameLhs
 
struct  OpElasticTools::OpRotatingFrameBoundaryLhs
 
struct  OpElasticTools::OpInternalForceRhs< LOGSTRAIN, REACTIONS >
 
struct  OpElasticTools::OpForceRhs
 
struct  OpElasticTools::OpCentrifugalForce
 
struct  OpElasticTools::OpStiffnessMatrixLhs
 
struct  OpElasticTools::OpLogStrainMatrixLhs
 
struct  OpElasticTools::OpPostProcElastic< LOGSTRAIN >
 
struct  OpElasticTools::OpSaveReactionForces
 [Operators definitions] More...
 
struct  OpElasticTools::OpBoundaryMassMatrix_UU
 
struct  OpElasticTools::OpEssentialRHS_U
 

Namespaces

 OpElasticTools
 

Macros

#define MAT_TO_TENSOR4_3D(MAT)
 

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)
 
auto get_rotation_R (Tensor1< double, 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 >
MoFEMErrorCode get_eigen_val_and_proj_souza (const Tensor2_symmetric< T, 3 > &X, double *e, Tensor2< double, 3, 3 > *E)
 
template<typename T >
auto get_ln_X (const Tensor2_symmetric< T, 3 > &X)
 
template<typename T >
auto calculate_d_05lnX_dX (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 >
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)
 
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)
 

Macro Definition Documentation

◆ MAT_TO_TENSOR4_3D

#define MAT_TO_TENSOR4_3D (   MAT)
Value:
&MAT(0, 0), &MAT(1, 0), &MAT(2, 0), &MAT(3, 0), &MAT(4, 0), &MAT(5, 0), \
&MAT(6, 0), &MAT(7, 0), &MAT(8, 0), &MAT(9, 0), &MAT(10, 0), \
&MAT(11, 0), &MAT(12, 0), &MAT(13, 0), &MAT(14, 0), &MAT(15, 0), \
&MAT(16, 0), &MAT(17, 0), &MAT(18, 0), &MAT(19, 0), &MAT(20, 0), \
&MAT(21, 0), &MAT(22, 0), &MAT(23, 0), &MAT(24, 0), &MAT(25, 0), \
&MAT(26, 0), &MAT(27, 0), &MAT(28, 0), &MAT(29, 0), &MAT(30, 0), \
&MAT(31, 0), &MAT(32, 0), &MAT(33, 0), &MAT(34, 0), &MAT(35, 0), \
&MAT(36, 0), &MAT(37, 0), &MAT(38, 0), &MAT(39, 0), &MAT(40, 0), \
&MAT(41, 0), &MAT(42, 0), &MAT(43, 0), &MAT(44, 0), &MAT(45, 0), \
&MAT(46, 0), &MAT(47, 0), &MAT(48, 0), &MAT(49, 0), &MAT(50, 0), \
&MAT(51, 0), &MAT(52, 0), &MAT(53, 0), &MAT(54, 0), &MAT(55, 0), \
&MAT(56, 0), &MAT(57, 0), &MAT(58, 0), &MAT(59, 0), &MAT(60, 0), \
&MAT(61, 0), &MAT(62, 0), &MAT(63, 0), &MAT(64, 0), &MAT(65, 0), \
&MAT(66, 0), &MAT(67, 0), &MAT(68, 0), &MAT(69, 0), &MAT(70, 0), \
&MAT(71, 0), &MAT(72, 0), &MAT(73, 0), &MAT(74, 0), &MAT(75, 0), \
&MAT(76, 0), &MAT(77, 0), &MAT(78, 0), &MAT(79, 0), &MAT(80, 0)

Definition at line 15 of file ElasticOperators.hpp.

Function Documentation

◆ calculate_d_05lnX_dX()

template<typename T >
auto calculate_d_05lnX_dX ( Tensor2_symmetric< T, 3 > &  X)

Definition at line 430 of file ElasticOperators.hpp.

430  {
431 
432  auto my_ln_x = get_diff_log_map(X);
433  my_ln_x(i, j, k, l) *= 0.5;
434  return my_ln_x;
435 
436  Tensor2<double, 3, 3> E[3];
437  double eig[3];
438  auto nX = to_non_symm(X);
439  CHKERR calcEigenvalues3x3(nX, eig);
440  CHKERR calcEigenproj3x3(nX, eig, E);
441  Tensor4<double, 3, 3, 3, 3> D_Xa;
442  CHKERR calcDerivativeLog3x3(nX, eig, E, D_Xa);
443  D_Xa(i, j, k, l) *= 0.5;
444  return tensor4_to_ddg(D_Xa);
445 };
MoFEMErrorCode calcEigenvalues3x3(const Tensor2< double, 3, 3 > &Tensor, double *eigen)
auto get_diff_log_map(Tensor2_symmetric< T, 3 > &X)
static Index< 'l', 3 > l
MoFEMErrorCode calcDerivativeLog3x3(const Tensor2< double, 3, 3 > &X, double *eigen, Tensor2< double, 3, 3 > *eigenProj, Tensor4< double, 3, 3, 3, 3 > &D_X)
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
static Index< 'i', 3 > i
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
static Index< 'k', 3 > k
auto to_non_symm(T &symm)
MoFEMErrorCode calcEigenproj3x3(const Tensor2< double, 3, 3 > &Tensor, double *eigen, Tensor2< double, 3, 3 > *eigenProj)

◆ calculate_nominal_moduli_TL()

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 
)

Definition at line 449 of file ElasticOperators.hpp.

453  {
454 
455  // Tensor4<double, 3, 3, 3, 3> TLs3;
456  Tensor4<double, 3, 3, 3, 3> P3;
457  Tensor2<double, 3, 3> invF;
458 
459  double det_f;
460  CHKERR determinantTensor3by3(F, det_f);
461  CHKERR invertTensor3by3(F, det_f, invF);
462 
463  // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
464 
465  // VectorDouble lambda(3);
466  VectorDouble e(3);
467  VectorDouble f(3);
468  VectorDouble d(3);
469  MatrixDouble Vi(3, 3);
470  MatrixDouble Ksi(3, 3);
471  MatrixDouble Zeta(3, 3);
472  Vi.clear();
473  Ksi.clear();
474  Zeta.clear();
475  double eta = 0;
476 
477  Tensor1<double, 3> _N[3];
478  Tensor1<double, 3> _n[3];
479 
480  for (int dd : {0, 1, 2}) {
481  for (int cc : {0, 1, 2})
482  _N[dd](cc) = eig_vec(dd, cc);
483  }
484 
485  for (int dd : {0, 1, 2})
486  _n[dd](i) = F(i, j) * _N[dd](j);
487 
488  for (int ii = 0; ii != 3; ++ii) {
489  e(ii) = 0.5 * log(lambda(ii));
490  d(ii) = 1. / lambda(ii);
491  f(ii) = -2. / (lambda(ii) * lambda(ii));
492  }
493 
494  const double &lam0 = lambda(0);
495  const double &lam1 = lambda(1);
496  const double &lam2 = lambda(2);
497 
498  if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
499  // all different
500  for (int ii = 0; ii != 3; ++ii)
501  for (int jj = 0; jj != 3; ++jj) {
502  if (ii != jj) {
503 
504  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
505  Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
506  for (int kk = 0; kk != 3; ++kk) {
507 
508  if (kk != ii && kk != jj)
509  eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
510  (lambda(ii) - lambda(kk)));
511  }
512  }
513  }
514 
515  } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
516  // all the same
517  for (int ii = 0; ii != 3; ++ii)
518  for (int jj = 0; jj != 3; ++jj)
519  if (ii != jj) {
520 
521  Vi(ii, jj) = 0.5 * d(ii);
522  Ksi(ii, jj) = 0.125 * f(ii);
523  }
524  eta = 0.125 * f(0);
525  } else if (is_eq(lam0, lam1)) {
526  int i_ = 0;
527  int j_ = 1;
528  int k_ = 2;
529  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
530  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
531 
532  for (int ii = 0; ii != 3; ++ii)
533  for (int jj = 0; jj != 3; ++jj)
534  if (ii != jj) {
535  if (jj == k_ || ii == k_) {
536  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
537  Ksi(ii, jj) =
538  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
539  }
540  }
541  eta = Ksi(k_, i_);
542  } else if (is_eq(lam0, lam2)) {
543  int i_ = 0;
544  int j_ = 2;
545  int k_ = 1;
546  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
547  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
548 
549  for (int ii = 0; ii != 3; ++ii)
550  for (int jj = 0; jj != 3; ++jj)
551  if (ii != jj) {
552  if (jj == k_ || ii == k_) {
553  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
554  Ksi(ii, jj) =
555  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
556  }
557  }
558  eta = Ksi(k_, i_);
559  } else if (is_eq(lam1, lam2)) {
560  int i_ = 1;
561  int j_ = 2;
562  int k_ = 0;
563  Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
564  Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
565 
566  for (int ii = 0; ii != 3; ++ii)
567  for (int jj = 0; jj != 3; ++jj)
568  if (ii != jj) {
569  if (jj == k_ || ii == k_) {
570  Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
571  Ksi(ii, jj) =
572  (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
573  }
574  }
575  eta = Ksi(k_, i_);
576  }
577 
578  P3(i, j, k, l) = 0.;
579 
580  Tensor2<double, 3, 3> Mii;
581  Tensor2<double, 3, 3> Mij;
582  Tensor2<double, 3, 3> Mik;
583  Tensor2<double, 3, 3> Mjk;
584  Tensor2<double, 3, 3> Mjj;
585 
586  // FIXME: TODO: put it nicely in one LOOP
587  for (int ii = 0; ii != 3; ++ii) {
588  Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
589  P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
590  for (int jj = 0; jj != 3; ++jj)
591  if (ii != jj) {
592  Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
593  P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
594  }
595  }
596 
597  if (!compute_tangent)
598  return P3;
599 
600  TLs3(i, j, k, l) = 0.;
601  for (int ii = 0; ii != 3; ++ii)
602  for (int jj = 0; jj != 3; ++jj)
603  Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
604 
605  for (int ii = 0; ii != 3; ++ii) {
606  Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
607  TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
608  }
609 
610  for (int ii = 0; ii != 3; ++ii)
611  for (int jj = 0; jj != 3; ++jj)
612  if (jj != ii)
613  for (int kk = 0; kk != 3; ++kk)
614  if (kk != ii && kk != jj) {
615  Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
616  Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
617  TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
618  }
619 
620  for (int ii = 0; ii != 3; ++ii)
621  for (int jj = 0; jj != 3; ++jj)
622  if (jj != ii) {
623  Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
624  Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
625  const double k2 = 2. * Ksi(ii, jj);
626  const double zp = k2 * Zeta(ii, jj);
627  const double zp2 = k2 * Zeta(jj, jj);
628  TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
629  zp * Mjj(i, j) * Mij(k, l) +
630  zp2 * Mij(i, j) * Mij(k, l);
631  }
632 
633  Tensor2<double, 3, 3> Piola;
634  Tensor2<double, 3, 3> invFPiola;
635  Tensor4<double, 3, 3, 3, 3> Ttmp3;
636  auto nstress = to_non_symm(stress);
637  Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
638  invFPiola(i, j) = invF(i, k) * Piola(k, j);
639  Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
640  0.5 * invFPiola(i, l) * kronecker_delta(j, k);
641  TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
642 
643  // TLs3(i, j, k, l) += invF(i, j) * (stress(m, n) * P3(m, n, k, l));
644 
645  // return TLs3;
646  return P3;
647 };
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:571
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
static Index< 'l', 3 > l
static Index< 'n', 3 > n
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:552
bool is_diff(const double &a, const double &b)
static Index< 'm', 3 > m
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static Index< 'i', 3 > i
bool is_eq(const double &a, const double &b)
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
static Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
auto to_non_symm(T &symm)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74

◆ exp_map()

template<typename T >
auto exp_map ( const Tensor2_symmetric< T, 3 > &  X)

Definition at line 154 of file ElasticOperators.hpp.

154  {
155  Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
156  Tensor2<double, 3, 3> X_nt;
157 
158  constexpr int max_it = 1000;
159  constexpr double eps = 1e-12;
160 
161  int fac_n = 1;
162 
163  auto Xcopy = to_non_symm(X);
164  auto X_n = to_non_symm(X);
165  exp_tensor(i, j) += Xcopy(i, j);
166 
167  for (int n = 2; n != max_it; ++n) {
168  fac_n *= n;
169  X_nt(i, j) = X_n(i, j);
170  X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
171  const double inv_nf = 1. / fac_n;
172  exp_tensor(i, j) += inv_nf * X_n(i, j);
173  const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
174  if (e < eps)
175  return to_symm(exp_tensor);
176  }
177 
178  return to_symm(exp_tensor);
179 };
auto to_symm(T &non_symm)
static Index< 'n', 3 > n
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k
auto to_non_symm(T &symm)
static const double eps

◆ get_diff_log_map()

template<typename T >
auto get_diff_log_map ( Tensor2_symmetric< T, 3 > &  X)

Definition at line 240 of file ElasticOperators.hpp.

240  {
241  Tensor2<double, 3, 3> X_nt_tmp;
242  Tensor4<double, 3, 3, 3, 3> diff_log_map;
243  Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
244 
245  constexpr int max_it = 1000;
246  constexpr double eps = 1e-12;
247 
248  auto get_Xn = [&](auto &log_tensor) {
249  vector<Tensor2<double, 3, 3>> Xn;
250  Xn.emplace_back();
251  auto &cu_Xn = Xn.back();
252  cu_Xn(i, j) = kronecker_delta(i, j);
253  for (int n = 2; n != max_it; ++n) {
254  const double cof = 1. / n;
255  auto X_nt_tmp = Xn.back();
256  Xn.emplace_back();
257  auto &c_Xn = Xn.back();
258  c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
259  const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
260  if (e < eps)
261  return Xn;
262  }
263  return Xn;
264  };
265 
266  auto log_tensor = to_non_symm(X);
267  log_tensor(i, j) -= kronecker_delta(i, j);
268  diff_log_map(i, j, k, l) = 0.;
269 
270  auto Xn_vec = get_Xn(log_tensor);
271 
272  for (int n = 1; n != Xn_vec.size() + 1; ++n) {
273  const double cof = pow(-1, n + 1) / n;
274  for (int m = 1; m != n + 1; ++m) {
275  auto &Xn_1 = Xn_vec[m - 1];
276  auto &Xn_2 = Xn_vec[n - m];
277  diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
278  }
279  }
280  // FIXME: make it more effient
281  auto isddg = tensor4_to_ddg((*cache).Is);
282  diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
283  return tensor4_to_ddg(diff_log_map_tmp);
284 };
static Index< 'l', 3 > l
static Index< 'n', 3 > n
static Index< 'm', 3 > m
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
auto to_non_symm(T &symm)
static const double eps

◆ get_eigen_val_and_proj_lapack()

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 
)

Definition at line 288 of file ElasticOperators.hpp.

290  {
291 
293  Tensor1<double, 3> eig;
294  Tensor2<double, 3, 3> eigen_vec_lap;
295  // eig_vec(i, j) = X(i, j);
296  eigen_vec_lap(0, 0) = X(0, 0);
297  eigen_vec_lap(1, 1) = X(1, 1);
298  eigen_vec_lap(2, 2) = X(2, 2);
299  eigen_vec_lap(0, 1) = eigen_vec_lap(1, 0) = X(0, 1);
300  eigen_vec_lap(2, 1) = eigen_vec_lap(1, 2) = X(2, 1);
301  eigen_vec_lap(2, 0) = eigen_vec_lap(0, 2) = X(2, 0);
302 
303  constexpr int dim = 3;
304  int n = dim, lda = dim, info, lwork = -1;
305 
306  double wkopt;
307  info = lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0), &wkopt,
308  lwork);
309  if (info != 0)
310  SETERRQ1(PETSC_COMM_SELF, 1,
311  "something is wrong with lapack_dsyev info = %d", info);
312 
313  lwork = (int)wkopt;
314  double work[lwork];
315  info = lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0), work,
316  lwork);
317  if (info != 0)
318  SETERRQ1(PETSC_COMM_SELF, 1,
319  "something is wrong with lapack_dsyev info = %d", info);
320  lambda(i) = eig(i);
321  eig_vec(i, j) = eigen_vec_lap(i, j);
322  // FIXME: perturb
323  // const double h = 1e-12;
324  // if (is_eq(lambda[0], lambda[1]))
325  // lambda[0] *= (1 + h);
326  // if (is_eq(lambda[1], lambda[2]))
327  // lambda[2] *= (1 - h);
328 
330 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'n', 3 > n
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
const int dim
static Index< 'i', 3 > i
static Index< 'j', 3 > j

◆ get_eigen_val_and_proj_souza()

template<typename T >
MoFEMErrorCode get_eigen_val_and_proj_souza ( const Tensor2_symmetric< T, 3 > &  X,
double e,
Tensor2< double, 3, 3 > *  E 
)

Definition at line 348 of file ElasticOperators.hpp.

349  {
350 
352  constexpr auto t_kd = Kronecker_Delta<int>();
353  constexpr double third = boost::math::constants::third<double>();
354  Tensor2<double, 3, 3> X2;
355  X2(i, j) = X(i, k) * X(k, j);
356 
357  const double I1 = X(i, i);
358  const double I2 = 0.5 * (I1 * I1 - X(j, i) * X(i, j));
359  double I3;
361 
362  const double R = (-2 * I1 * I1 * I1 + 9 * I1 * I2 - 27 * I3) / 54.;
363  const double Q = (I1 * I1 - 3 * I2) / 9.;
364  const double theta = acos(R / sqrt(Q * Q * Q));
365  double sqrtQ = sqrt(Q);
366  if (isnormal(sqrtQ)) {
367  e[0] = -2 * sqrtQ * cos(theta * third) + I1 * third;
368  e[1] = -2 * sqrtQ * cos((theta + 2 * M_PI) * third) + I1 * third;
369  e[2] = -2 * sqrtQ * cos((theta - 2 * M_PI) * third) + I1 * third;
370  } else {
371  e[0] = I1 * third;
372  e[1] = I1 * third;
373  e[2] = I1 * third;
374  }
375 
376  if (is_diff(e[0], e[1]) && is_diff(e[1], e[2]) && is_diff(e[0], e[2])) {
377  for (int ii = 0; ii != 3; ++ii) {
378  const double &x = e[ii];
379  const double coef = x / (2 * x * x * x - I1 * x * x + I3);
380  E[ii](i, j) = coef * (X2(i, j) - (I1 - x) * X(i, j) +
381  (I3 / x) * kronecker_delta(i, j));
382  }
383 
384  } else if (is_eq(e[0], e[1]) && is_eq(e[1], e[2])) {
385  E[0](i, j) = t_kd(i, j);
386  E[1](i, j) = 0.;
387  E[2](i, j) = 0.;
388  } else if (is_eq(e[1], e[2])) {
389  E[0](i, j) =
390  (e[0] / (2. * e[0] * e[0] * e[0] - I1 * e[0] * e[0] + I3)) *
391  (X(i, k) * X(k, j) - (I1 - e[0]) * X(i, j) + (I3 / e[0]) * t_kd(i, j));
392  E[1](i, j) = t_kd(i, j) - E[0](i, j);
393  E[2](i, j) = 0.;
394  } else if (is_eq(e[0], e[1])) {
395  E[2](i, j) =
396  (e[2] / (2. * e[2] * e[2] * e[2] - I1 * e[2] * e[2] + I3)) *
397  (X(i, k) * X(k, j) - (I1 - e[2]) * X(i, j) + (I3 / e[2]) * t_kd(i, j));
398  E[0](i, j) = t_kd(i, j) - E[2](i, j);
399  E[1](i, j) = 0.;
400  } else {
401  E[1](i, j) =
402  (e[1] / (2. * e[1] * e[1] * e[1] - I1 * e[1] * e[1] + I3)) *
403  (X(i, k) * X(k, j) - (I1 - e[1]) * X(i, j) + (I3 / e[1]) * t_kd(i, j));
404  E[2](i, j) = t_kd(i, j) - E[1](i, j);
405  E[0](i, j) = 0.;
406  }
407 
409 }
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:571
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
bool is_diff(const double &a, const double &b)
static Index< 'i', 3 > i
bool is_eq(const double &a, const double &b)
static Index< 'j', 3 > j
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double R
static Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
constexpr double third

◆ get_ln_X()

template<typename T >
auto get_ln_X ( const Tensor2_symmetric< T, 3 > &  X)

Definition at line 411 of file ElasticOperators.hpp.

411  {
412 
413  // Tensor2<double, 3, 3> lnX(0., 0., 0., 0., 0., 0., 0., 0., 0.);
414  // Tensor2<double, 3, 3> E[3];
415  // double lambda[3];
416  // auto nX = to_non_symm(X);
417  // CHKERR calcEigenvalues3x3(nX, lambda);
418  // CHKERR calcEigenproj3x3(nX, lambda, E);
419  // // CHKERR get_eigen_val_and_proj_souza(X, lambda, E);
420 
421  // for (int ii = 0; ii != 3; ii++)
422  // lnX(i, j) += E[ii](i, j) * log(lambda[ii]);
423 
424  auto test = get_log_map(X);
425  return get_log_map(X);
426  // return to_symm(lnX);
427 };
auto get_log_map(const Tensor2_symmetric< T, 3 > &X)

◆ get_ln_X_lapack()

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 
)

Definition at line 333 of file ElasticOperators.hpp.

335  {
336  Tensor2_symmetric<double, 3> lnX;
337  Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
338  for (int ii = 0; ii != 3; ++ii)
339  diag(ii, ii) = log(lambda(ii));
340 
341  lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
342 
343  return lnX;
344 };
static Index< 'l', 3 > l
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ get_log_map()

template<typename T >
auto get_log_map ( const Tensor2_symmetric< T, 3 > &  X)

Definition at line 214 of file ElasticOperators.hpp.

214  {
215 
216  constexpr int max_it = 1000;
217  constexpr double eps = 1e-12;
218 
219  auto log_tensor = to_non_symm(X);
220  log_tensor(i, j) -= kronecker_delta(i, j);
221 
222  auto Xcopy = to_non_symm(log_tensor);
223  auto X_n = to_non_symm(log_tensor);
224  Tensor2<double, 3, 3> X_nt_tmp;
225 
226  for (int n = 2; n != max_it; ++n) {
227 
228  X_nt_tmp(i, j) = X_n(i, j);
229  X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
230  const double inv_n = 1. / n;
231  log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
232  const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
233  if (e < eps)
234  break;
235  }
236 
237  return to_symm(log_tensor);
238 };
auto to_symm(T &non_symm)
static Index< 'n', 3 > n
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
auto to_non_symm(T &symm)
static const double eps

◆ get_rotation_R()

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

Definition at line 181 of file ElasticOperators.hpp.

181  {
182 
184 
185  auto get_rotation = [&](FTensor1 &t_omega) {
187 
188  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
189  t_R(i, j) = t_kd(i, j);
190 
191  const double angle = sqrt(t_omega(i) * t_omega(i));
192  if (std::abs(angle) < 1e-18)
193  return t_R;
194 
196  t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
197  const double a = sin(angle) / angle;
198  const double ss_2 = sin(angle / 2.);
199  const double b = 2. * ss_2 * ss_2 / (angle * angle);
200  t_R(i, j) += a * t_Omega(i, j);
201  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
202 
203  return t_R;
204  };
205 
206  Tensor1<double, 3> my_c_omega;
207  my_c_omega(i) = t1_omega(i) * tt;
208  auto t_R = get_rotation(my_c_omega);
209 
210  return get_rotation(my_c_omega);
211 };
Kronecker Delta class.
static Index< 'i', 3 > i
static Index< 'j', 3 > j
static Index< 'k', 3 > k

◆ get_tensor4_from_ddg()

template<typename T >
auto get_tensor4_from_ddg ( Ddg< T, 3, 3 > &  ddg)

Definition at line 129 of file ElasticOperators.hpp.

129  {
130  Tensor4<double, 3, 3, 3, 3> tens4;
131  for (int ii = 0; ii != 3; ++ii)
132  for (int jj = 0; jj != 3; ++jj)
133  for (int kk = 0; kk != 3; ++kk)
134  for (int ll : {0, 1, 2})
135  tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
136  return tens4;
137 };

◆ is_ddg()

template<typename T >
bool is_ddg ( Tensor4< T, 3, 3, 3, 3 >  tensor)

Definition at line 139 of file ElasticOperators.hpp.

139  {
140  auto test_ddg = tensor4_to_ddg(tensor);
141  for (int ii = 0; ii != 3; ++ii)
142  for (int jj = 0; jj != 3; ++jj)
143  for (int kk = 0; kk != 3; ++kk)
144  for (int ll : {0, 1, 2}) {
145  const double diff =
146  abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
147  if (diff > 1e-6)
148  return false;
149  }
150 
151  return true;
152 };
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)

◆ is_diff()

bool is_diff ( const double a,
const double b 
)

Definition at line 76 of file ElasticOperators.hpp.

76  {
77  constexpr double eps = std::numeric_limits<double>::epsilon();
78  return abs(a - b) > eps;
79 };
static const double eps

◆ is_eq()

bool is_eq ( const double a,
const double b 
)

Definition at line 72 of file ElasticOperators.hpp.

72  {
73  constexpr double eps = std::numeric_limits<double>::epsilon();
74  return abs(a - b) < eps;
75 };
static const double eps

◆ tensor4_to_ddg()

template<typename T >
Ddg<double, 3, 3> tensor4_to_ddg ( Tensor4< T, 3, 3, 3, 3 > &  t)

Definition at line 82 of file ElasticOperators.hpp.

82  {
83  Ddg<double, 3, 3> loc_tens;
84 
85  Number<0> N0;
86  Number<1> N1;
87  Number<2> N2;
88 
89  loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
90  loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
91  loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
92  loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
93  loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
94  loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
95  loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
96  loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
97  loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
98  loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
99  loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
100  loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
101  loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
102  loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
103  loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
104  loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
105  loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
106  loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
107  loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
108  loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
109  loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
110  loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
111  loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
112  loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
113  loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
114  loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
115  loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
116  loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
117  loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
118  loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
119  loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
120  loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
121  loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
122  loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
123  loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
124  loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
125 
126  return loc_tens;
127 };
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0

◆ to_non_symm()

template<typename T >
auto to_non_symm ( T symm)

Definition at line 33 of file ElasticOperators.hpp.

33  {
34  Tensor2<double, 3, 3> non_symm;
35  Number<0> N0;
36  Number<1> N1;
37  Number<2> N2;
38  non_symm(N0, N0) = symm(N0, N0);
39  non_symm(N1, N1) = symm(N1, N1);
40  non_symm(N2, N2) = symm(N2, N2);
41  non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
42  non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
43  non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
44  return non_symm;
45 };
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0

◆ to_symm()

template<typename T >
auto to_symm ( T non_symm)

Definition at line 47 of file ElasticOperators.hpp.

47  {
48  Tensor2_symmetric<double, 3> symm;
49  Number<0> N0;
50  Number<1> N1;
51  Number<2> N2;
52  symm(N0, N0) = non_symm(N0, N0);
53  symm(N1, N1) = non_symm(N1, N1);
54  symm(N2, N2) = non_symm(N2, N2);
55  symm(N0, N1) = non_symm(N0, N1);
56  symm(N0, N2) = non_symm(N0, N2);
57  symm(N1, N2) = non_symm(N1, N2);
58  return symm;
59 };
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0