v0.13.2
Loading...
Searching...
No Matches
Classes | Typedefs | Functions
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< OpMassPrecHTensorBoundOpSchurPreconditionMassContactBound
 
typedef struct OpSchurPreconditionMassContact< OpMassPrecHTensorVolOpSchurPreconditionMassContactVol
 

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

Definition at line 1194 of file ElasticOperators.hpp.

◆ 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.

613 {
614
615 Ddg<double, 3, 3> P3;
616 // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
617
618 VectorDouble e(3);
619 VectorDouble f(3);
620 VectorDouble d(3);
621 MatrixDouble Vi(3, 3);
622 MatrixDouble Ksi(3, 3);
623 MatrixDouble Zeta(3, 3);
624 Vi.clear();
625 Ksi.clear();
626 Zeta.clear();
627 double eta = 0;
628
629 Tensor1<double, 3> _N[3];
630 for (int dd = 0; dd != 3; ++dd) {
631 for (int cc = 0; cc != 3; ++cc)
632 _N[dd](cc) = eig_vec(dd, cc);
633 }
634
635 for (int ii = 0; ii != 3; ++ii) {
636 e(ii) = 0.5 * log(lambda(ii));
637 d(ii) = 1. / lambda(ii);
638 f(ii) = -2. / (lambda(ii) * lambda(ii));
639 }
640
641 const double &lam0 = lambda(0);
642 const double &lam1 = lambda(1);
643 const double &lam2 = lambda(2);
644
645 auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
646 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
647 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
648
649 for (int ii = 0; ii != 3; ++ii)
650 for (int jj = 0; jj != 3; ++jj)
651 if (ii != jj) {
652 if (jj == k_ || ii == k_) {
653 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
654 Ksi(ii, jj) =
655 (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
656 }
657 }
658 eta = Ksi(k_, i_);
659 };
660
661 if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
662 // all different
663 for (int ii = 0; ii != 3; ++ii)
664 for (int jj = 0; jj != 3; ++jj) {
665 if (ii != jj) {
666
667 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
668 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
669 for (int kk = 0; kk != 3; ++kk) {
670
671 if (kk != ii && kk != jj)
672 eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
673 (lambda(ii) - lambda(kk)));
674 }
675 }
676 }
677
678 } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
679 // all the same
680 for (int ii = 0; ii != 3; ++ii)
681 for (int jj = 0; jj != 3; ++jj)
682 if (ii != jj) {
683
684 Vi(ii, jj) = 0.5 * d(ii);
685 Ksi(ii, jj) = 0.125 * f(ii);
686 }
687 eta = 0.125 * f(0);
688 } else if (is_eq(lam0, lam1)) {
689
690 compute_ksi_eta_2eq(0, 1, 2);
691
692 } else if (is_eq(lam0, lam2)) {
693
694 compute_ksi_eta_2eq(0, 2, 1);
695
696 } else if (is_eq(lam1, lam2)) {
697
698 compute_ksi_eta_2eq(1, 2, 0);
699 }
700
701 P3(i, j, k, l) = 0.;
702 Tensor2_symmetric<double, 3> Mii;
703 Tensor2_symmetric<double, 3> Mij;
704 Tensor2_symmetric<double, 3> Mik;
705 Tensor2_symmetric<double, 3> Mjk;
706 Tensor2_symmetric<double, 3> Mjj;
707
708 for (int ii = 0; ii != 3; ++ii) {
709 Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
710 P3(i, j, k, l) += ((0.5 * d(ii) * _N[ii](i)) ^ _N[ii](j)) * Mii(k, l);
711 for (int jj = 0; jj != 3; ++jj)
712 if (ii != jj) {
713 Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
714 P3(i, j, k, l) += ((Vi(ii, jj) * _N[ii](i)) ^ _N[jj](j)) * Mij(k, l);
715 }
716 }
717
718 if (!TANGENT)
719 return P3;
720
721 TLs3(i, j, k, l) = 0.;
722 for (int ii = 0; ii != 3; ++ii)
723 for (int jj = 0; jj != 3; ++jj)
724 Zeta(ii, jj) = stress(i, j) * (_N[ii](i) ^ _N[jj](j));
725
726 for (int ii = 0; ii != 3; ++ii) {
727 Mii(i, j) = (_N[ii](i) ^ _N[ii](j)) + (_N[ii](i) ^ _N[ii](j));
728 TLs3(i, j, k, l) += (0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j)) * Mii(k, l);
729 }
730
731 for (int ii = 0; ii != 3; ++ii)
732 for (int jj = 0; jj != 3; ++jj)
733 if (jj != ii) {
734 Mij(i, j) = (_N[ii](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[ii](j));
735 Mjj(i, j) = (_N[jj](i) ^ _N[jj](j)) + (_N[jj](i) ^ _N[jj](j));
736 const double k2 = 2. * Ksi(ii, jj);
737 const double zp = k2 * Zeta(ii, jj);
738 const double zp2 = k2 * Zeta(jj, jj);
739 TLs3(i, j, k, l) += (zp * Mij(i, j)) * Mjj(k, l) +
740 (zp * Mjj(i, j)) * Mij(k, l) +
741 (zp2 * Mij(i, j)) * Mij(k, l);
742 for (int kk = 0; kk != 3; ++kk) {
743 if (kk != ii && kk != jj) {
744 Mik(i, j) = (_N[ii](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[ii](j));
745 Mjk(i, j) = (_N[jj](i) ^ _N[kk](j)) + (_N[kk](i) ^ _N[jj](j));
746 TLs3(i, j, k, l) +=
747 (2. * eta * Zeta(ii, jj) * Mik(i, j)) * Mjk(k, l);
748 }
749 }
750 }
751
752 return P3;
753};
constexpr double lambda
surface tension
constexpr double eta
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
bool is_eq(const double &a, const double &b)
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
bool is_diff(const double &a, const double &b)

◆ 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.

386 {
387
388 // Tensor4<double, 3, 3, 3, 3> TLs3;
389 Tensor4<double, 3, 3, 3, 3> P3;
390 Tensor2<double, 3, 3> invF;
391
392 double det_f;
393 CHKERR determinantTensor3by3(F, det_f);
394 CHKERR invertTensor3by3(F, det_f, invF);
395
396 // from https://www.sciencedirect.com/science/article/pii/S0045782502004383
397
398 // VectorDouble lambda(3);
399 VectorDouble e(3);
400 VectorDouble f(3);
401 VectorDouble d(3);
402 MatrixDouble Vi(3, 3);
403 MatrixDouble Ksi(3, 3);
404 MatrixDouble Zeta(3, 3);
405 Vi.clear();
406 Ksi.clear();
407 Zeta.clear();
408 double eta = 0;
409
410 Tensor1<double, 3> _N[3];
411 Tensor1<double, 3> _n[3];
412
413 for (int dd : {0, 1, 2}) {
414 for (int cc : {0, 1, 2})
415 _N[dd](cc) = eig_vec(dd, cc);
416 }
417
418 for (int dd : {0, 1, 2})
419 _n[dd](i) = F(i, j) * _N[dd](j);
420
421 for (int ii = 0; ii != 3; ++ii) {
422 e(ii) = 0.5 * log(lambda(ii));
423 d(ii) = 1. / lambda(ii);
424 f(ii) = -2. / (lambda(ii) * lambda(ii));
425 }
426
427 const double &lam0 = lambda(0);
428 const double &lam1 = lambda(1);
429 const double &lam2 = lambda(2);
430
431 auto compute_ksi_eta_2eq = [&](int i_, int j_, int k_) {
432 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
433 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
434
435 for (int ii = 0; ii != 3; ++ii)
436 for (int jj = 0; jj != 3; ++jj)
437 if (ii != jj) {
438 if (jj == k_ || ii == k_) {
439 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
440 Ksi(ii, jj) =
441 (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
442 }
443 }
444 eta = Ksi(k_, i_);
445 };
446
447 if (is_diff(lam0, lam1) && is_diff(lam0, lam2) && is_diff(lam1, lam2)) {
448 // all different
449 for (int ii = 0; ii != 3; ++ii)
450 for (int jj = 0; jj != 3; ++jj) {
451 if (ii != jj) {
452
453 Vi(ii, jj) = (e(ii) - e(jj)) / (lambda(ii) - lambda(jj));
454 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (lambda(ii) - lambda(jj));
455 for (int kk = 0; kk != 3; ++kk) {
456
457 if (kk != ii && kk != jj)
458 eta += e(ii) / (2. * (lambda(ii) - lambda(jj)) *
459 (lambda(ii) - lambda(kk)));
460 }
461 }
462 }
463
464 } else if (is_eq(lam0, lam1) && is_eq(lam0, lam2) && is_eq(lam1, lam2)) {
465 // all the same
466 for (int ii = 0; ii != 3; ++ii)
467 for (int jj = 0; jj != 3; ++jj)
468 if (ii != jj) {
469
470 Vi(ii, jj) = 0.5 * d(ii);
471 Ksi(ii, jj) = 0.125 * f(ii);
472 }
473 eta = 0.125 * f(0);
474 } else if (is_eq(lam0, lam1)) {
475
476 compute_ksi_eta_2eq(0, 1, 2);
477
478 } else if (is_eq(lam0, lam2)) {
479
480 compute_ksi_eta_2eq(0, 2, 1);
481
482 } else if (is_eq(lam1, lam2)) {
483
484 compute_ksi_eta_2eq(1, 2, 0);
485 }
486
487 P3(i, j, k, l) = 0.;
488
489 Tensor2<double, 3, 3> Mii;
490 Tensor2<double, 3, 3> Mij;
491 Tensor2<double, 3, 3> Mik;
492 Tensor2<double, 3, 3> Mjk;
493 Tensor2<double, 3, 3> Mjj;
494
495 for (int ii = 0; ii != 3; ++ii) {
496 Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
497 P3(i, j, k, l) += 0.5 * d(ii) * _N[ii](i) * _N[ii](j) * Mii(k, l);
498 for (int jj = 0; jj != 3; ++jj)
499 if (ii != jj) {
500 Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
501 P3(i, j, k, l) += Vi(ii, jj) * _N[ii](i) * _N[jj](j) * Mij(k, l);
502 }
503 }
504
505 if (!compute_tangent)
506 return P3;
507
508 TLs3(i, j, k, l) = 0.;
509 for (int ii = 0; ii != 3; ++ii)
510 for (int jj = 0; jj != 3; ++jj)
511 Zeta(ii, jj) = (stress(i, j) * (_N[ii](i) * _N[jj](j)));
512
513 for (int ii = 0; ii != 3; ++ii) {
514 Mii(i, j) = _n[ii](i) * _N[ii](j) + _n[ii](i) * _N[ii](j);
515 TLs3(i, j, k, l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(i, j) * Mii(k, l);
516 }
517
518 for (int ii = 0; ii != 3; ++ii)
519 for (int jj = 0; jj != 3; ++jj)
520 if (jj != ii)
521 for (int kk = 0; kk != 3; ++kk)
522 if (kk != ii && kk != jj) {
523 Mik(i, j) = _n[ii](i) * _N[kk](j) + _n[kk](i) * _N[ii](j);
524 Mjk(i, j) = _n[jj](i) * _N[kk](j) + _n[kk](i) * _N[jj](j);
525 TLs3(i, j, k, l) += 2. * eta * Zeta(ii, jj) * Mik(i, j) * Mjk(k, l);
526 }
527
528 for (int ii = 0; ii != 3; ++ii)
529 for (int jj = 0; jj != 3; ++jj)
530 if (jj != ii) {
531 Mij(i, j) = _n[ii](i) * _N[jj](j) + _n[jj](i) * _N[ii](j);
532 Mjj(i, j) = _n[jj](i) * _N[jj](j) + _n[jj](i) * _N[jj](j);
533 const double k2 = 2. * Ksi(ii, jj);
534 const double zp = k2 * Zeta(ii, jj);
535 const double zp2 = k2 * Zeta(jj, jj);
536 TLs3(i, j, k, l) += zp * Mij(i, j) * Mjj(k, l) +
537 zp * Mjj(i, j) * Mij(k, l) +
538 zp2 * Mij(i, j) * Mij(k, l);
539 }
540
541 Tensor2<double, 3, 3> Piola;
542 Tensor2<double, 3, 3> invFPiola;
543 Tensor4<double, 3, 3, 3, 3> Ttmp3;
544 auto nstress = to_non_symm(stress);
545 Piola(k, l) = nstress(m, n) * P3(m, n, k, l);
546 invFPiola(i, j) = invF(i, k) * Piola(k, j);
547 Ttmp3(i, j, k, l) = 0.5 * invFPiola(i, k) * kronecker_delta(j, l) +
548 0.5 * invFPiola(i, l) * kronecker_delta(j, k);
549 TLs3(i, j, k, l) += Ttmp3(i, j, k, l);
550
551 return P3;
552};
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
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
auto to_non_symm(T &symm)

◆ exp_map()

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

Definition at line 200 of file ElasticOperators.hpp.

200 {
201 Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
202 Tensor2<double, 3, 3> X_nt;
203
204 constexpr int max_it = 1000;
205 constexpr double eps = 1e-12;
206
207 int fac_n = 1;
208
209 auto Xcopy = to_non_symm(X);
210 auto X_n = to_non_symm(X);
211 exp_tensor(i, j) += Xcopy(i, j);
212
213 for (int n = 2; n != max_it; ++n) {
214 fac_n *= n;
215 X_nt(i, j) = X_n(i, j);
216 X_n(i, k) = X_nt(i, j) * Xcopy(j, k);
217 const double inv_nf = 1. / fac_n;
218 exp_tensor(i, j) += inv_nf * X_n(i, j);
219 const double e = inv_nf * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
220 if (e < eps)
221 return to_symm(exp_tensor);
222 }
223
224 return to_symm(exp_tensor);
225};
static const double eps
auto to_symm(T &non_symm)

◆ 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.

285 {
286 Tensor2<double, 3, 3> X_nt_tmp;
287 Tensor4<double, 3, 3, 3, 3> diff_log_map;
288 Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
289
290 constexpr int max_it = 1000;
291 constexpr double eps = 1e-12;
292
293 auto get_Xn = [&](auto &log_tensor) {
294 vector<Tensor2<double, 3, 3>> Xn;
295 Xn.emplace_back();
296 auto &cu_Xn = Xn.back();
297 cu_Xn(i, j) = kronecker_delta(i, j);
298 for (int n = 2; n != max_it; ++n) {
299 const double cof = 1. / n;
300 auto X_nt_tmp = Xn.back();
301 Xn.emplace_back();
302 auto &c_Xn = Xn.back();
303 c_Xn(i, k) = X_nt_tmp(i, j) * log_tensor(j, k);
304 const double e = cof * std::abs(sqrt(c_Xn(i, j) * c_Xn(i, j)));
305 if (e < eps)
306 return Xn;
307 }
308 return Xn;
309 };
310
311 auto log_tensor = to_non_symm(X);
312 log_tensor(i, j) -= kronecker_delta(i, j);
313 diff_log_map(i, j, k, l) = 0.;
314
315 auto Xn_vec = get_Xn(log_tensor);
316
317 for (int n = 1; n != Xn_vec.size() + 1; ++n) {
318 const double cof = pow(-1, n + 1) / n;
319 for (int m = 1; m != n + 1; ++m) {
320 auto &Xn_1 = Xn_vec[m - 1];
321 auto &Xn_2 = Xn_vec[n - m];
322 diff_log_map(i, j, k, l) += cof * Xn_1(i, k) * Xn_2(l, j);
323 }
324 }
325
326 auto isddg = tensor4_to_ddg((*cache).Is);
327 diff_log_map_tmp(i, j, k, l) = isddg(i, j, m, n) * diff_log_map(m, n, k, l);
328 return tensor4_to_ddg(diff_log_map_tmp);
329};
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)

◆ 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.

335 {
336
339 auto eigen_vec_lap = to_non_symm(X);
340 constexpr int n = 3;
341 constexpr int lda = 3;
342 constexpr int lwork = 15;
343 std::array<double, 15> work;
344 if (lapack_dsyev('V', 'U', n, &(eigen_vec_lap(0, 0)), lda, &eig(0),
345 work.data(), lwork) > 0)
346 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
347 "The algorithm failed to compute eigenvalues.");
348
349 lambda(i) = eig(i);
350 eig_vec(i, j) = eigen_vec_lap(i, j);
351
352 // perturb
353 // const double h = 1e-12;
354 // if (is_eq(lambda[0], lambda[1]))
355 // lambda[0] *= (1 + h);
356 // if (is_eq(lambda[1], lambda[2]))
357 // lambda[2] *= (1 - h);
358
360}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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:261

◆ 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.

376 {
377 return get_log_map(X);
378};
auto get_log_map(const Tensor2_symmetric< T, 3 > &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.

365 {
366 Tensor2_symmetric<double, 3> lnX;
367 Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
368 for (int ii = 0; ii != 3; ++ii)
369 diag(ii, ii) = log(lambda(ii));
370
371 lnX(i, j) = eig_vec(l, i) ^ (diag(l, k) * eig_vec(k, j));
372
373 return lnX;
374};

◆ 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.

259 {
260
261 constexpr int max_it = 1000;
262 constexpr double eps = 1e-12;
263
264 auto log_tensor = to_non_symm(X);
265 log_tensor(i, j) -= kronecker_delta(i, j);
266
267 auto Xcopy = to_non_symm(log_tensor);
268 auto X_n = to_non_symm(log_tensor);
269 Tensor2<double, 3, 3> X_nt_tmp;
270
271 for (int n = 2; n != max_it; ++n) {
272
273 X_nt_tmp(i, j) = X_n(i, j);
274 X_n(i, k) = X_nt_tmp(i, j) * Xcopy(j, k);
275 const double inv_n = 1. / n;
276 log_tensor(i, j) += inv_n * pow(-1, n + 1) * X_n(i, j);
277 const double e = inv_n * std::abs(sqrt(X_n(i, j) * X_n(i, j)));
278 if (e < eps)
279 break;
280 }
281
282 return to_symm(log_tensor);
283};

◆ 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.

228 {
229
230 auto get_rotation = [&](Tensor1<double, 3> &t_omega) {
232
233 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
234 t_R(i, j) = t_kd(i, j);
235
236 const double angle = sqrt(t_omega(i) * t_omega(i));
237 if (std::abs(angle) < 1e-18)
238 return t_R;
239
241 t_Omega(i, j) = FTensor::levi_civita<double>(i, j, k) * t_omega(k);
242 const double a = sin(angle) / angle;
243 const double ss_2 = sin(angle / 2.);
244 const double b = 2. * ss_2 * ss_2 / (angle * angle);
245 t_R(i, j) += a * t_Omega(i, j);
246 t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
247
248 return t_R;
249 };
250
251 Tensor1<double, 3> my_c_omega;
252 my_c_omega(i) = t1_omega(i) * tt;
253 auto t_R = get_rotation(my_c_omega);
254
255 return get_rotation(my_c_omega);
256};
constexpr double a
Kronecker Delta class.
constexpr auto t_kd

◆ 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.

175 {
176 Tensor4<double, 3, 3, 3, 3> tens4;
177 for (int ii = 0; ii != 3; ++ii)
178 for (int jj = 0; jj != 3; ++jj)
179 for (int kk = 0; kk != 3; ++kk)
180 for (int ll = 0; ll != 3; ++ll)
181 tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
182 return tens4;
183};

◆ 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.

818 {
819 Tensor4<double, 3, 3, 3, 3> my_D;
820 for (int k = 0; k != 3; ++k) {
821 for (int l = 0; l != 3; ++l) {
822 auto d_stress = getFiniteDiff(k, l, common_data, F, t_stress);
823 for (int i = 0; i != 3; ++i) {
824 for (int j = 0; j != 3; ++j) {
825 my_D(i, j, k, l) = d_stress(i, j);
826 }
827 }
828 }
829 }
830 return my_D;
831};
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)

◆ 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.

758 {
760 Tensor2<double, 3, 3> out;
761
762 Tensor2<double, 3, 3> strainPl;
763 Tensor2<double, 3, 3> strainMin;
764 strainPl(i, j) = F(i, j);
765 strainMin(i, j) = F(i, j);
766 const double h = 1e-8;
767 strainPl(kk, ll) += h;
768 strainMin(kk, ll) -= h;
769 Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
770
771 auto calculate_stress = [&](auto &F) {
772 Tensor2<double, 3, 3> Piola;
773 Tensor2_symmetric<double, 3> stress_sym;
774 Tensor2_symmetric<double, 3> C;
775 Tensor2_symmetric<double, 3> strain;
776
777 C(i, j) = F(m, i) ^ F(m, j);
778
780 Tensor2<double, 3, 3> eig_vec;
781 CHKERR get_eigen_val_and_proj_lapack(C, lambda, eig_vec);
782
783 auto lnC = get_ln_X_lapack(C, lambda, eig_vec);
784 // auto lnC = get_log_map(C);
785 strain(i, j) = 0.5 * lnC(i, j);
786 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*common_data.mtD);
787 stress_sym(i, j) = t_D(i, j, k, l) * strain(k, l);
788 auto stress_tmp = to_non_symm(stress_sym);
789
790 // Tensor4<double, 3, 3, 3, 3> D2;
791 // constexpr auto t_kd = Kronecker_Delta<int>();
792 // Tensor4<double, 3, 3, 3, 3> dC_dF;
793 // dC_dF(i, j, k, l) = (t_kd(i, l) * F(k, j)) + (t_kd(j, l) * F(k, i));
794 // auto t_L = get_diff_log_map(C);
795 // t_L(i, j, k, l) *= 0.5;
796 // D2(i, j, k, l) = t_L(i, j, m, n) * dC_dF(m, n, k, l);
797
798 auto D2 = calculate_nominal_moduli_TL(C, F, stress_sym, lambda, eig_vec,
799 dummy, false);
800
801 Piola(k, l) = stress_tmp(i, j) * D2(i, j, k, l);
802
803 return Piola;
804 };
805
806 auto StressPlusH = calculate_stress(strainPl);
807 auto StressMinusH = calculate_stress(strainMin);
808 // out(i, j) = (StressPlusH(i, j) - t_stress(i, j)) / (h);
809
810 out(i, j) = (StressPlusH(i, j) - StressMinusH(i, j)) / (2 * h);
811
812 return out;
813}
auto get_ln_X_lapack(const Tensor2_symmetric< T1, 3 > &X, Tensor1< T2, 3 > &lambda, Tensor2< T3, 3, 3 > &eig_vec)
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)
double h
boost::shared_ptr< MatrixDouble > mtD

◆ 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.

185 {
186 auto test_ddg = tensor4_to_ddg(tensor);
187 for (int ii = 0; ii != 3; ++ii)
188 for (int jj = 0; jj != 3; ++jj)
189 for (int kk = 0; kk != 3; ++kk)
190 for (int ll : {0, 1, 2}) {
191 const double diff =
192 abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
193 if (diff > 1e-6)
194 return false;
195 }
196
197 return true;
198};

◆ is_diff()

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

Definition at line 122 of file ElasticOperators.hpp.

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

◆ is_eq()

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

Definition at line 118 of file ElasticOperators.hpp.

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

◆ 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.

128 {
129 Ddg<double, 3, 3> loc_tens;
130
131 Number<0> N0;
132 Number<1> N1;
133 Number<2> N2;
134
135 loc_tens(N0, N0, N0, N0) = t(N0, N0, N0, N0); // 0
136 loc_tens(N0, N0, N0, N1) = t(N0, N0, N0, N1); // 1
137 loc_tens(N0, N0, N0, N2) = t(N0, N0, N0, N2); // 2
138 loc_tens(N0, N0, N1, N1) = t(N0, N0, N1, N1); // 3
139 loc_tens(N0, N0, N1, N2) = t(N0, N0, N1, N2); // 4
140 loc_tens(N0, N0, N2, N2) = t(N0, N0, N2, N2); // 5
141 loc_tens(N0, N1, N0, N0) = t(N0, N1, N0, N0); // 6
142 loc_tens(N0, N1, N0, N1) = t(N0, N1, N0, N1); // 7
143 loc_tens(N0, N1, N0, N2) = t(N0, N1, N0, N2); // 8
144 loc_tens(N0, N1, N1, N1) = t(N0, N1, N1, N1); // 9
145 loc_tens(N0, N1, N1, N2) = t(N0, N1, N1, N2); // 10
146 loc_tens(N0, N1, N2, N2) = t(N0, N1, N2, N2); // 11
147 loc_tens(N0, N2, N0, N0) = t(N0, N2, N0, N0); // 12
148 loc_tens(N0, N2, N0, N1) = t(N0, N2, N0, N1); // 13
149 loc_tens(N0, N2, N0, N2) = t(N0, N2, N0, N2); // 14
150 loc_tens(N0, N2, N1, N1) = t(N0, N2, N1, N1); // 15
151 loc_tens(N0, N2, N1, N2) = t(N0, N2, N1, N2); // 16
152 loc_tens(N0, N2, N2, N2) = t(N0, N2, N2, N2); // 17
153 loc_tens(N1, N1, N0, N0) = t(N1, N1, N0, N0); // 18
154 loc_tens(N1, N1, N0, N1) = t(N1, N1, N0, N1); // 19
155 loc_tens(N1, N1, N0, N2) = t(N1, N1, N0, N2); // 20
156 loc_tens(N1, N1, N1, N1) = t(N1, N1, N1, N1); // 21
157 loc_tens(N1, N1, N1, N2) = t(N1, N1, N1, N2); // 22
158 loc_tens(N1, N1, N2, N2) = t(N1, N1, N2, N2); // 23
159 loc_tens(N1, N2, N0, N0) = t(N1, N2, N0, N0); // 24
160 loc_tens(N1, N2, N0, N1) = t(N1, N2, N0, N1); // 25
161 loc_tens(N1, N2, N0, N2) = t(N1, N2, N0, N2); // 26
162 loc_tens(N1, N2, N1, N1) = t(N1, N2, N1, N1); // 27
163 loc_tens(N1, N2, N1, N2) = t(N1, N2, N1, N2); // 28
164 loc_tens(N1, N2, N2, N2) = t(N1, N2, N2, N2); // 29
165 loc_tens(N2, N2, N0, N0) = t(N2, N2, N0, N0); // 30
166 loc_tens(N2, N2, N0, N1) = t(N2, N2, N0, N1); // 31
167 loc_tens(N2, N2, N0, N2) = t(N2, N2, N0, N2); // 32
168 loc_tens(N2, N2, N1, N1) = t(N2, N2, N1, N1); // 33
169 loc_tens(N2, N2, N1, N2) = t(N2, N2, N1, N2); // 34
170 loc_tens(N2, N2, N2, N2) = t(N2, N2, N2, N2); // 35
171
172 return loc_tens;
173};
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
constexpr double t
plate stiffness
Definition: plate.cpp:59

◆ 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.

558 {
559
560 Tensor1<double, 3> eig{lambda(0), lambda(1), lambda(2)};
561 Tensor2<double, 3, 3> eigen_vec;
562 eigen_vec(i, j) = eig_vec(i, j);
563
564 auto get_uniq_nb = [&](auto eig) {
565 return distance(&eig(0), unique(&eig(0), &eig(0) + 3, is_eq));
566 };
567
568 auto sort_eigen_vals2 = [&](auto &eig, auto &eigen_vec) {
569 if (is_eq(eig(0), eig(1))) {
570 Tensor2<double, 3, 3> eigen_vec_c{
571 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
572 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
573 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
574 Tensor1<double, 3> eig_c{eig(0), eig(2), eig(1)};
575 eig(i) = eig_c(i);
576 eigen_vec(i, j) = eigen_vec_c(i, j);
577 } else {
578 Tensor2<double, 3, 3> eigen_vec_c{
579 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
580 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
581 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
582 Tensor1<double, 3> eig_c{eig(1), eig(0), eig(2)};
583 eig(i) = eig_c(i);
584 eigen_vec(i, j) = eigen_vec_c(i, j);
585 }
586 };
587
588 Tensor2_symmetric<double, 3> Ts;
589 Ts(i, j) = stress(i, j);
590 auto f = [](double v) { return 0.5 * log(v); };
591 auto d_f = [](double v) { return 0.5 / v; };
592 auto dd_f = [](double v) { return -0.5 / (v * v); };
593
594 auto nb_uniq = get_uniq_nb(eig);
595 if (nb_uniq == 2)
596 sort_eigen_vals2(eig, eigen_vec);
597 auto t_d =
598
599 EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
600
601 auto t_dd =
602
603 EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, Ts, nb_uniq);
604
605 TLs3(i, j, k, l) = t_dd(i, j, k, l);
606 return t_d;
607};
const double v
phase velocity of light in medium (cm/ns)
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto dd_f
Definition: HenckyOps.hpp:17
auto d_f
Definition: HenckyOps.hpp:16

◆ to_non_symm()

template<typename T >
auto OpElasticTools::to_non_symm ( T symm)
inline

Definition at line 90 of file ElasticOperators.hpp.

90 {
91 Tensor2<double, 3, 3> non_symm;
92 Number<0> N0;
93 Number<1> N1;
94 Number<2> N2;
95 non_symm(N0, N0) = symm(N0, N0);
96 non_symm(N1, N1) = symm(N1, N1);
97 non_symm(N2, N2) = symm(N2, N2);
98 non_symm(N0, N1) = non_symm(N1, N0) = symm(N0, N1);
99 non_symm(N0, N2) = non_symm(N2, N0) = symm(N0, N2);
100 non_symm(N2, N1) = non_symm(N1, N2) = symm(N2, N1);
101 return non_symm;
102};

◆ to_symm()

template<typename T >
auto OpElasticTools::to_symm ( T non_symm)
inline

Definition at line 104 of file ElasticOperators.hpp.

104 {
105 Tensor2_symmetric<double, 3> symm;
106 Number<0> N0;
107 Number<1> N1;
108 Number<2> N2;
109 symm(N0, N0) = non_symm(N0, N0);
110 symm(N1, N1) = non_symm(N1, N1);
111 symm(N2, N2) = non_symm(N2, N2);
112 symm(N0, N1) = non_symm(N0, N1);
113 symm(N0, N2) = non_symm(N0, N2);
114 symm(N1, N2) = non_symm(N1, N2);
115 return symm;
116};