493                                              {
  495 
  502  
  504 
  505  auto nb_gauss_pts = DomainEleOp::getGaussPts().size2();
  506  auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
  508  auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
 
  509  auto t_f = getFTensor0FromVec(
commonDataPtr->plasticSurface);
 
  510  auto t_flow = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticFlow);
 
  511  auto t_plastic_strain =
  512      getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticStrain);
 
  513  auto t_plastic_strain_dot =
  514      getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->plasticStrainDot);
 
  515  auto t_stress =
  516      getFTensor2SymmetricFromMat<DIM>(*(
commonDataPtr->mStressPtr));
 
  517 
  518  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
commonDataPtr->mDPtr);
 
  519  auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*
mDPtr);
 
  520 
  523 
  526  t_flow_dir_dstress(
i, 
j, 
k, 
l) =
 
  527      1.5 * (t_diff_deviator(
M, 
N, 
i, 
j) * t_diff_deviator(
M, 
N, 
k, 
l));
 
  528  t_flow_dir_dstrain(
i, 
j, 
k, 
l) =
 
  529      t_flow_dir_dstress(
i, 
j, 
m, 
n) * t_D_Op(
m, 
n, 
k, 
l);
 
  530 
  531 
  532  auto t_alpha_dir =
  534 
  542                                       false);
  544                                          false);
  545 
  554 
  556  auto t_res_c_dtau = getFTensor0FromVec(
commonDataPtr->resCdTau);
 
  557  auto t_res_c_dstrain =
  558      getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resCdStrain);
 
  559  auto t_res_c_plastic_strain =
  560      getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resCdPlasticStrain);
 
  561  auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resFlow);
 
  562  auto t_res_flow_dtau =
  563      getFTensor2SymmetricFromMat<DIM>(
commonDataPtr->resFlowDtau);
 
  564  auto t_res_flow_dstrain =
  565      getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->resFlowDstrain);
 
  566  auto t_res_flow_dplastic_strain =
  567      getFTensor4DdgFromMat<DIM, DIM>(
commonDataPtr->resFlowDstrainDot);
 
  568 
  569  auto next = [&]() {
  570    ++t_tau;
  571    ++t_tau_dot;
  572    ++t_f;
  573    ++t_flow;
  574    ++t_plastic_strain;
  575    ++t_plastic_strain_dot;
  576    ++t_stress;
  577    ++t_res_c;
  578    ++t_res_c_dtau;
  579    ++t_res_c_dstrain;
  580    ++t_res_c_plastic_strain;
  581    ++t_res_flow;
  582    ++t_res_flow_dtau;
  583    ++t_res_flow_dstrain;
  584    ++t_res_flow_dplastic_strain;
  585    ++t_w;
  586  };
  587 
  588  auto get_avtive_pts = [&]() {
  589    int nb_points_avtive_on_elem = 0;
  590    int nb_points_on_elem = 0;
  591 
  593    auto t_tau_dot = getFTensor0FromVec(
commonDataPtr->plasticTauDot);
 
  594    auto t_f = getFTensor0FromVec(
commonDataPtr->plasticSurface);
 
  595    auto t_plastic_strain_dot =
  596        getFTensor2SymmetricFromMat<SPACE_DIM>(
commonDataPtr->plasticStrainDot);
 
  597 
  598    auto dt = this->getTStimeStep();
 
  599 
  600    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
  603          eqiv, t_tau_dot, t_f,
  608 
  609      ++nb_points_on_elem;
  610      if (sign_ww > 0) {
  611        ++nb_points_avtive_on_elem;
  612      }
  613 
  614      ++t_tau;
  615      ++t_tau_dot;
  616      ++t_f;
  617      ++t_plastic_strain_dot;
  618    }
  619 
  625 
  626    ++nb_elements;
  627    nb_points += nb_points_on_elem;
  628    if (nb_points_avtive_on_elem > 0) {
  629      ++avtive_elems;
  630      active_points += nb_points_avtive_on_elem;
  631      if (nb_points_avtive_on_elem == nb_points_on_elem) {
  632        ++avtive_full_elems;
  633      }
  634    }
  635 
  636    if (nb_points_avtive_on_elem != nb_points_on_elem)
  637      return 1;
  638    else
  639      return 0;
  640  };
  641 
  642  if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
  643    get_avtive_pts();
  644  }
  645 
  646  auto dt = this->getTStimeStep();
 
  647  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
  648 
  651                                                  t_diff_plastic_strain,
  653 
  654    const auto sigma_y =
  657    const auto d_sigma_y =
  660 
  664 
  665    auto c = 
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
 
  674 
  676 
  677        t_stress, 
trace(t_stress),
 
  678 
  680 
  681    );
  682 
  684    t_flow_dir(
k, 
l) = 1.5 * (t_dev_stress(
I, 
J) * t_diff_deviator(
I, 
J, 
k, 
l));
 
  686    t_flow_dstrain(
i, 
j) = t_flow(
k, 
l) * t_D_Op(
k, 
l, 
i, 
j);
 
  687 
  688    auto get_res_c = [&]() { 
return c; };
 
  689 
  690    auto get_res_c_dstrain = [&](auto &t_diff_res) {
  691      t_diff_res(
i, 
j) = c_f * t_flow_dstrain(
i, 
j);
 
  692    };
  693 
  694    auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
  695      t_diff_res(
i, 
j) = (this->getTSa() * c_equiv) * t_diff_eqiv(
i, 
j);
 
  696      t_diff_res(
k, 
l) -= c_f * t_flow(
i, 
j) * t_alpha_dir(
i, 
j, 
k, 
l);
 
  697    };
  698 
  699    auto get_res_c_dtau = [&]() {
  700      return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
  701    };
  702 
  703    auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
  704      t_diff_res(
k, 
l) = -c_f * t_flow(
i, 
j) * t_alpha_dir(
i, 
j, 
k, 
l);
 
  705    };
  706 
  707    auto get_res_flow = [&](auto &t_res_flow) {
  708      const auto a = sigma_y;
 
  709      const auto b = t_tau_dot;
  710      t_res_flow(
k, 
l) = 
a * t_plastic_strain_dot(
k, 
l) - b * t_flow_dir(
k, 
l);
 
  711    };
  712 
  713    auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
  714      const auto da = d_sigma_y;
  715      const auto db = this->getTSa();
  716      t_res_flow_dtau(
k, 
l) =
 
  717          da * t_plastic_strain_dot(
k, 
l) - db * t_flow_dir(
k, 
l);
 
  718    };
  719 
  720    auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
  721      const auto b = t_tau_dot;
  722      t_res_flow_dstrain(
m, 
n, 
k, 
l) = -t_flow_dir_dstrain(
m, 
n, 
k, 
l) * b;
 
  723    };
  724 
  725    auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
  726      const auto a = sigma_y;
 
  727      t_res_flow_dplastic_strain(
m, 
n, 
k, 
l) =
 
  728          (
a * this->getTSa()) * t_diff_plastic_strain(
m, 
n, 
k, 
l);
 
  729      const auto b = t_tau_dot;
  730      t_res_flow_dplastic_strain(
m, 
n, 
i, 
j) +=
 
  731          (t_flow_dir_dstrain(
m, 
n, 
k, 
l) * t_alpha_dir(
k, 
l, 
i, 
j)) * b;
 
  732    };
  733 
  734    t_res_c = get_res_c();
  735    get_res_flow(t_res_flow);
  736 
  737    if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
  738      t_res_c_dtau = get_res_c_dtau();
  739      get_res_c_dstrain(t_res_c_dstrain);
  740      get_res_c_dplastic_strain(t_res_c_plastic_strain);
  741      get_res_flow_dtau(t_res_flow_dtau);
  742      get_res_flow_dstrain(t_res_flow_dstrain);
  743      get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
  744    }
  745 
  746    next();
  747  }
  748 
  750}
#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
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< 'M', 3 > M
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
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)
FTensor::Index< 'N', 3 > N
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
[Common data]
FTensor::Index< 'J', 3 > J
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
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)
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)
FTensor::Index< 'm', 3 > m
static std::array< int, 5 > activityData