492 {
494
501
503
504 auto nb_gauss_pts = DomainEleOp::getGaussPts().size2();
505 auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
507 auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
508 auto t_f = getFTensor0FromVec(
commonDataPtr->plasticSurface);
509 auto t_flow = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticFlow);
510 auto t_plastic_strain =
511 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticStrain);
512 auto t_plastic_strain_dot =
513 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticStrainDot);
514 auto t_stress = getFTensor2SymmetricFromMat<DIM>(*(
commonDataPtr->mStressPtr));
515
516 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
commonDataPtr->mDPtr);
517 auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*
mDPtr);
518
521
524 t_flow_dir_dstress(
i,
j,
k,
l) =
525 1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l));
526 t_flow_dir_dstrain(
i,
j,
k,
l) =
527 t_flow_dir_dstress(
i,
j,
m,
n) * t_D_Op(
m,
n,
k,
l);
528
529
530 auto t_alpha_dir =
532
540 false);
542 false);
543
552
554 auto t_res_c_dtau = getFTensor0FromVec(
commonDataPtr->resCdTau);
555 auto t_res_c_dstrain =
556 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resCdStrain);
557 auto t_res_c_plastic_strain =
558 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resCdPlasticStrain);
559 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resFlow);
560 auto t_res_flow_dtau =
561 getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resFlowDtau);
562 auto t_res_flow_dstrain =
563 getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->resFlowDstrain);
564 auto t_res_flow_dplastic_strain =
565 getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->resFlowDstrainDot);
566
567 auto next = [&]() {
568 ++t_tau;
569 ++t_tau_dot;
570 ++t_f;
571 ++t_flow;
572 ++t_plastic_strain;
573 ++t_plastic_strain_dot;
574 ++t_stress;
575 ++t_res_c;
576 ++t_res_c_dtau;
577 ++t_res_c_dstrain;
578 ++t_res_c_plastic_strain;
579 ++t_res_flow;
580 ++t_res_flow_dtau;
581 ++t_res_flow_dstrain;
582 ++t_res_flow_dplastic_strain;
583 ++t_w;
584 };
585
586 auto get_avtive_pts = [&]() {
587 int nb_points_avtive_on_elem = 0;
588 int nb_points_on_elem = 0;
589
591 auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
592 auto t_f = getFTensor0FromVec(
commonDataPtr->plasticSurface);
593 auto t_plastic_strain_dot =
594 getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
595
596 auto dt = this->getTStimeStep();
597
598 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
601 eqiv, t_tau_dot, t_f,
606
607 ++nb_points_on_elem;
608 if (sign_ww > 0) {
609 ++nb_points_avtive_on_elem;
610 }
611
612 ++t_tau;
613 ++t_tau_dot;
614 ++t_f;
615 ++t_plastic_strain_dot;
616 }
617
623
624 ++nb_elements;
625 nb_points += nb_points_on_elem;
626 if (nb_points_avtive_on_elem > 0) {
627 ++avtive_elems;
628 active_points += nb_points_avtive_on_elem;
629 if (nb_points_avtive_on_elem == nb_points_on_elem) {
630 ++avtive_full_elems;
631 }
632 }
633
634 if (nb_points_avtive_on_elem != nb_points_on_elem)
635 return 1;
636 else
637 return 0;
638 };
639
640 if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
641 get_avtive_pts();
642 }
643
644 auto dt = this->getTStimeStep();
645 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
646
649 t_diff_plastic_strain,
651
652 const auto sigma_y =
655 const auto d_sigma_y =
658
662
663 auto c =
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
672
674
675 t_stress,
trace(t_stress),
676
678
679 );
680
682 t_flow_dir(
k,
l) = 1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l));
684 t_flow_dstrain(
i,
j) = t_flow(
k,
l) * t_D_Op(
k,
l,
i,
j);
685
686 auto get_res_c = [&]() {
return c; };
687
688 auto get_res_c_dstrain = [&](auto &t_diff_res) {
689 t_diff_res(
i,
j) = c_f * t_flow_dstrain(
i,
j);
690 };
691
692 auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
693 t_diff_res(
i,
j) = (this->getTSa() * c_equiv) * t_diff_eqiv(
i,
j);
694 t_diff_res(
k,
l) -= c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
695 };
696
697 auto get_res_c_dtau = [&]() {
698 return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
699 };
700
701 auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
702 t_diff_res(
k,
l) = -c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
703 };
704
705 auto get_res_flow = [&](auto &t_res_flow) {
706 const auto a = sigma_y;
707 const auto b = t_tau_dot;
708 t_res_flow(
k,
l) =
a * t_plastic_strain_dot(
k,
l) - b * t_flow_dir(
k,
l);
709 };
710
711 auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
712 const auto da = d_sigma_y;
713 const auto db = this->getTSa();
714 t_res_flow_dtau(
k,
l) =
715 da * t_plastic_strain_dot(
k,
l) - db * t_flow_dir(
k,
l);
716 };
717
718 auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
719 const auto b = t_tau_dot;
720 t_res_flow_dstrain(
m,
n,
k,
l) = -t_flow_dir_dstrain(
m,
n,
k,
l) * b;
721 };
722
723 auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
724 const auto a = sigma_y;
725 t_res_flow_dplastic_strain(
m,
n,
k,
l) =
726 (
a * this->getTSa()) * t_diff_plastic_strain(
m,
n,
k,
l);
727 const auto b = t_tau_dot;
728 t_res_flow_dplastic_strain(
m,
n,
i,
j) +=
729 (t_flow_dir_dstrain(
m,
n,
k,
l) * t_alpha_dir(
k,
l,
i,
j)) * b;
730 };
731
732 t_res_c = get_res_c();
733 get_res_flow(t_res_flow);
734
735 if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
736 t_res_c_dtau = get_res_c_dtau();
737 get_res_c_dstrain(t_res_c_dstrain);
738 get_res_c_dplastic_strain(t_res_c_plastic_strain);
739 get_res_flow_dtau(t_res_flow_dtau);
740 get_res_flow_dstrain(t_res_flow_dstrain);
741 get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
742 }
743
744 next();
745 }
746
748}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Index< 'N', 3 > N
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
double constrian_sign(double x, double dt)
FTensor::Index< 'J', 3 > J
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
FTensor::Index< 'M', 3 > M
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain, FTensor::Number< DIM >)
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w, double vis_H, double sigma_Y)
auto diff_constrain_df(double sign)
auto diff_constrain_dsigma_y(double sign)
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
FTensor::Index< 'I', 3 > I
[Common data]
auto diff_deviator(FTensor::Ddg< double, DIM, DIM > &&t_diff_stress, FTensor::Number< DIM >)
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
auto deviator(FTensor::Tensor2_symmetric< T, DIM > &t_stress, double trace, FTensor::Tensor2_symmetric< double, DIM > &t_alpha, FTensor::Number< DIM >)
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
double constrain_abs(double x, double dt)
FTensor::Index< 'm', 3 > m
static std::array< int, 5 > activityData
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
