52  t_diff(
i, 
j, 
k, 
l) = 0;
 
   53  t_diff(0, 0, 0, 0) = 1;
 
   54  t_diff(1, 1, 1, 1) = 1;
 
   56  t_diff(1, 0, 1, 0) = 0.5;
 
   57  t_diff(1, 0, 0, 1) = 0.5;
 
   59  t_diff(0, 1, 0, 1) = 0.5;
 
   60  t_diff(0, 1, 1, 0) = 0.5;
 
   62  if constexpr (DIM == 3) {
 
   63    t_diff(2, 2, 2, 2) = 1;
 
   65    t_diff(2, 0, 2, 0) = 0.5;
 
   66    t_diff(2, 0, 0, 2) = 0.5;
 
   67    t_diff(0, 2, 0, 2) = 0.5;
 
   68    t_diff(0, 2, 2, 0) = 0.5;
 
   70    t_diff(2, 1, 2, 1) = 0.5;
 
   71    t_diff(2, 1, 1, 2) = 0.5;
 
   72    t_diff(1, 2, 1, 2) = 0.5;
 
   73    t_diff(1, 2, 2, 1) = 0.5;
 
 
  130  t_diff_deviator(
I, 
J, 
k, 
l) = 0;
 
  131  for (
int ii = 0; ii != DIM; ++ii)
 
  132    for (
int jj = ii; jj != DIM; ++jj)
 
  133      for (
int kk = 0; kk != DIM; ++kk)
 
  134        for (
int ll = kk; ll != DIM; ++ll)
 
  135          t_diff_deviator(ii, jj, kk, ll) = t_diff_stress(ii, jj, kk, ll);
 
  137  constexpr double third = boost::math::constants::third<double>();
 
  139  t_diff_deviator(0, 0, 0, 0) -= 
third;
 
  140  t_diff_deviator(0, 0, 1, 1) -= 
third;
 
  142  t_diff_deviator(1, 1, 0, 0) -= 
third;
 
  143  t_diff_deviator(1, 1, 1, 1) -= 
third;
 
  145  t_diff_deviator(2, 2, 0, 0) -= 
third;
 
  146  t_diff_deviator(2, 2, 1, 1) -= 
third;
 
  148  if constexpr (DIM == 3) {
 
  149    t_diff_deviator(0, 0, 2, 2) -= 
third;
 
  150    t_diff_deviator(1, 1, 2, 2) -= 
third;
 
  151    t_diff_deviator(2, 2, 2, 2) -= 
third;
 
  154  return t_diff_deviator;
 
 
  388      &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
 
  389      &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
 
  390      &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
 
  391      &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
 
  392      &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
 
  393      &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
 
 
  493    int side, EntityType type, 
EntData &data) {
 
  503  auto ¶ms = commonDataPtr->blockParams; 
 
  505  auto nb_gauss_pts = DomainEleOp::getGaussPts().size2();
 
  506  auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
 
  507  auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
 
  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);
 
  516      getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
 
  518  auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
 
  519  auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
 
  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);
 
  535  commonDataPtr->resC.resize(nb_gauss_pts, 
false);
 
  536  commonDataPtr->resCdTau.resize(nb_gauss_pts, 
false);
 
  537  commonDataPtr->resCdStrain.resize(
size_symm, nb_gauss_pts, 
false);
 
  538  commonDataPtr->resCdPlasticStrain.resize(
size_symm, nb_gauss_pts, 
false);
 
  539  commonDataPtr->resFlow.resize(
size_symm, nb_gauss_pts, 
false);
 
  540  commonDataPtr->resFlowDtau.resize(
size_symm, nb_gauss_pts, 
false);
 
  546  commonDataPtr->resC.clear();
 
  547  commonDataPtr->resCdTau.clear();
 
  548  commonDataPtr->resCdStrain.clear();
 
  549  commonDataPtr->resCdPlasticStrain.clear();
 
  550  commonDataPtr->resFlow.clear();
 
  551  commonDataPtr->resFlowDtau.clear();
 
  552  commonDataPtr->resFlowDstrain.clear();
 
  553  commonDataPtr->resFlowDstrainDot.clear();
 
  555  auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
 
  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);
 
  575    ++t_plastic_strain_dot;
 
  580    ++t_res_c_plastic_strain;
 
  583    ++t_res_flow_dstrain;
 
  584    ++t_res_flow_dplastic_strain;
 
  588  auto get_avtive_pts = [&]() {
 
  589    int nb_points_avtive_on_elem = 0;
 
  590    int nb_points_on_elem = 0;
 
  592    auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
 
  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);
 
  598    auto dt = this->getTStimeStep();
 
  600    for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
  603          eqiv, t_tau_dot, t_f,
 
  611        ++nb_points_avtive_on_elem;
 
  617      ++t_plastic_strain_dot;
 
  627    nb_points += nb_points_on_elem;
 
  628    if (nb_points_avtive_on_elem > 0) {
 
  630      active_points += nb_points_avtive_on_elem;
 
  631      if (nb_points_avtive_on_elem == nb_points_on_elem) {
 
  636    if (nb_points_avtive_on_elem != nb_points_on_elem)
 
  642  if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
 
  646  auto dt = this->getTStimeStep();
 
  647  for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
  651                                                  t_diff_plastic_strain,
 
  657    const auto d_sigma_y =
 
  665    auto c = 
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
 
  677        t_stress, 
trace(t_stress),
 
  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);
 
  688    auto get_res_c = [&]() { 
return c; };
 
  690    auto get_res_c_dstrain = [&](
auto &t_diff_res) {
 
  691      t_diff_res(
i, 
j) = c_f * t_flow_dstrain(
i, 
j);
 
  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);
 
  699    auto get_res_c_dtau = [&]() {
 
  700      return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
 
  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);
 
  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);
 
  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);
 
  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;
 
  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;
 
  734    t_res_c = get_res_c();
 
  735    get_res_flow(t_res_flow);
 
  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);
 
 
  984      &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
 
  985      &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
 
  986      &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
 
  987      &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
 
  988      &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
 
  989      &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
 
  990      &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
 
  991      &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
 
  992      &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
 
  993      &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
 
  994      &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
 
  995      &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
 
 
 1001    EntitiesFieldData::EntData &row_data,
 
 1002    EntitiesFieldData::EntData &col_data) {
 
 1009  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
 
 1013  auto &locMat = AssemblyDomainEleOp::locMat;
 
 1015  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
 
 1016  const auto nb_row_base_functions = row_data.getN().size2();
 
 1018  auto t_res_flow_dstrain =
 
 1019      getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
 
 1020  auto t_res_flow_dplastic_strain =
 
 1021      getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
 
 1025    ++t_res_flow_dstrain;
 
 1026    ++t_res_flow_dplastic_strain;
 
 1029  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
 
 1030  auto t_row_base = row_data.getFTensor0N();
 
 1031  for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
 1032    double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
 
 1037        alpha * (t_L(
i, 
j, O) * ((t_res_flow_dplastic_strain(
i, 
j, 
k, 
l) -
 
 1038                                  t_res_flow_dstrain(
i, 
j, 
k, 
l)) *
 
 1043    for (; rr != AssemblyDomainEleOp::nbRows / 
size_symm; ++rr) {
 
 1046      auto t_col_base = col_data.getFTensor0N(gg, 0);
 
 1047      for (
size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / 
size_symm; ++cc) {
 
 1048        t_mat(O, L) += ((t_row_base * t_col_base) * t_res_mat(O, L));
 
 1056    for (; rr < nb_row_base_functions; ++rr)
 
 
 1121    EntitiesFieldData::EntData &row_data,
 
 1122    EntitiesFieldData::EntData &col_data) {
 
 1127  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
 
 1130  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
 
 1131  const size_t nb_row_base_functions = row_data.getN().size2();
 
 1132  auto &locMat = AssemblyDomainEleOp::locMat;
 
 1134  auto t_res_flow_dtau =
 
 1135      getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
 
 1139  auto next = [&]() { ++t_res_flow_dtau; };
 
 1141  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
 
 1142  auto t_row_base = row_data.getFTensor0N();
 
 1143  for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
 1144    double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
 
 1147    t_res_vec(L) = alpha * (t_res_flow_dtau(
i, 
j) * t_L(
i, 
j, L));
 
 1151    for (; rr != AssemblyDomainEleOp::nbRows / 
size_symm; ++rr) {
 
 1154      auto t_col_base = col_data.getFTensor0N(gg, 0);
 
 1155      for (
size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; cc++) {
 
 1156        t_mat(L) += t_row_base * t_col_base * t_res_vec(L);
 
 1162    for (; rr != nb_row_base_functions; ++rr)
 
 
 1209    EntitiesFieldData::EntData &row_data,
 
 1210    EntitiesFieldData::EntData &col_data) {
 
 1215  constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
 
 1218  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
 
 1219  const auto nb_row_base_functions = row_data.getN().size2();
 
 1222      getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
 
 1223  auto t_c_dplastic_strain =
 
 1224      getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
 
 1228    ++t_c_dplastic_strain;
 
 1233  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
 
 1234  auto t_row_base = row_data.getFTensor0N();
 
 1235  for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
 
 1236    const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
 
 1241        t_L(
i, 
j, L) * (t_c_dplastic_strain(
i, 
j) - t_c_dstrain(
i, 
j));
 
 1247    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
 
 1248      const auto row_base = alpha * t_row_base;
 
 1249      auto t_col_base = col_data.getFTensor0N(gg, 0);
 
 1250      for (
size_t cc = 0; cc != AssemblyDomainEleOp::nbCols / 
size_symm; cc++) {
 
 1251        t_mat(L) += (row_base * t_col_base) * t_res_vec(L);
 
 1257    for (; rr != nb_row_base_functions; ++rr)
 
 
 1291    EntitiesFieldData::EntData &row_data,
 
 1292    EntitiesFieldData::EntData &col_data) {
 
 1295  const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
 
 1296  const auto nb_row_base_functions = row_data.getN().size2();
 
 1298  auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
 
 1299  auto next = [&]() { ++t_res_c_dtau; };
 
 1301  auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
 
 1302  auto t_row_base = row_data.getFTensor0N();
 
 1303  for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
 1304    const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
 
 1307    const auto res = alpha * (t_res_c_dtau);
 
 1310    auto mat_ptr = AssemblyDomainEleOp::locMat.data().begin();
 
 1312    for (; rr != AssemblyDomainEleOp::nbRows; ++rr) {
 
 1313      auto t_col_base = col_data.getFTensor0N(gg, 0);
 
 1314      for (
size_t cc = 0; cc != AssemblyDomainEleOp::nbCols; ++cc) {
 
 1315        *mat_ptr += t_row_base * t_col_base * res;
 
 1321    for (; rr < nb_row_base_functions; ++rr)