24 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
28 if constexpr (DIM == 2) {
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;
81 constexpr double third = boost::math::constants::third<double>();
82 return (t_stress(0, 0) + t_stress(1, 1)) *
third;
87 constexpr double third = boost::math::constants::third<double>();
88 return (t_stress(0, 0) + t_stress(1, 1) + t_stress(2, 2)) *
third;
91template <
typename T,
int DIM>
100 for (
int ii = 0; ii != DIM; ++ii)
101 for (
int jj = ii; jj != DIM; ++jj)
102 t_dev(ii, jj) = t_stress(ii, jj);
103 t_dev(0, 0) -=
trace;
104 t_dev(1, 1) -=
trace;
105 t_dev(2, 2) -=
trace;
106 for (
int ii = 0; ii != DIM; ++ii)
107 for (
int jj = ii; jj != DIM; ++jj)
108 t_dev(ii, jj) -= t_alpha(ii, jj);
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;
190 return std::sqrt(1.5 * t_stress_deviator(
I,
J) * t_stress_deviator(
I,
J)) +
191 std::numeric_limits<double>::epsilon();
202 (1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l))) / f;
206template <
typename T,
int DIM>
216 t_diff_flow(
i,
j,
k,
l) =
217 (1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l) -
218 (2. / 3.) * t_flow(
i,
j) * t_flow(
k,
l))) /
223template <
typename T,
int DIM>
234 t_diff_flow(
i,
j,
k,
l) =
235 t_diff_plastic_flow_dstress(
i,
j,
m,
n) * t_D(
m,
n,
k,
l);
252 const auto y = x /
zeta;
253 if (y > std::numeric_limits<float>::max_exponent10 ||
254 y < std::numeric_limits<float>::min_exponent10) {
260 const auto e = std::exp(y);
261 return (e - 1) / (1 + e);
266 const auto y = -x /
zeta;
267 if (y > std::numeric_limits<float>::max_exponent10 ||
268 y < std::numeric_limits<float>::min_exponent10) {
271 const double e = std::exp(y);
272 return x + 2 *
zeta * std::log1p(e);
276inline double w(
double eqiv,
double dot_tau,
double f,
double sigma_y,
278 return (f - sigma_y) / sigma_Y +
cn1 * (eqiv);
292inline double constraint(
double eqiv,
double dot_tau,
double f,
double sigma_y,
293 double abs_w,
double vis_H,
double sigma_Y) {
294 return vis_H * dot_tau +
295 (sigma_Y / 2) * ((
cn0 * (dot_tau - eqiv) +
cn1 * (eqiv) -
296 (f - sigma_y) / sigma_Y) -
301 double vis_H,
double sigma_Y) {
302 return vis_H + (sigma_Y / 2) * (
cn0);
307 return (sigma_Y / 2) * (-
cn0 +
cn1 * (1 - sign));
314template <
typename T,
int DIM>
322 return t_diff_constrain_dstress;
325template <
typename T1,
typename T2,
int DIM>
332 t_diff_constrain_dstrain(
k,
l) =
333 t_diff_constrain_dstress(
i,
j) * t_D(
i,
j,
k,
l);
334 return t_diff_constrain_dstrain;
337template <
typename T,
int DIM>
342 constexpr double A = 2. / 3;
343 return std::sqrt(
A * t_plastic_strain_dot(
i,
j) *
344 t_plastic_strain_dot(
i,
j)) +
345 std::numeric_limits<double>::epsilon();
348template <
typename T1,
typename T2,
typename T3,
int DIM>
350 T3 &t_diff_plastic_strain,
356 constexpr double A = 2. / 3;
358 t_diff_eqiv(
i,
j) =
A * (t_plastic_strain_dot(
k,
l) / eqiv) *
359 t_diff_plastic_strain(
k,
l,
i,
j);
369 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
370 &mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
376 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
377 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
378 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
379 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
380 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
381 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2)};
386template <
int DIM,
typename DomainEleOp>
390 boost::shared_ptr<CommonData> common_data_ptr);
397template <
int DIM,
typename DomainEleOp>
400 boost::shared_ptr<CommonData> common_data_ptr)
402 commonDataPtr(common_data_ptr) {
408template <
int DIM,
typename DomainEleOp>
416 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
418 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
419 auto t_plastic_strain =
420 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
422 commonDataPtr->plasticSurface.resize(nb_gauss_pts,
false);
423 commonDataPtr->plasticFlow.resize(
size_symm, nb_gauss_pts,
false);
424 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
426 auto ¶ms = commonDataPtr->blockParams;
428 for (
auto &f : commonDataPtr->plasticSurface) {
433 t_stress,
trace(t_stress),
446 t_flow(
i,
j) = t_flow_tmp(
i,
j);
456template <
int DIM,
typename DomainEleOp>
459 boost::shared_ptr<CommonData> common_data_ptr,
460 boost::shared_ptr<MatrixDouble> m_D_ptr);
465 boost::shared_ptr<MatrixDouble>
mDPtr;
468template <
int DIM,
typename DomainEleOp>
470 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
471 boost::shared_ptr<MatrixDouble> m_D_ptr)
473 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
479template <
int DIM,
typename DomainEleOp>
491 auto ¶ms = commonDataPtr->blockParams;
495 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
496 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
497 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
498 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
499 auto t_plastic_strain =
500 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
501 auto t_plastic_strain_dot =
502 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
504 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
506 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
507 auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
514 t_flow_dir_dstress(
i,
j,
k,
l) =
515 1.5 * (t_diff_deviator(
M,
N,
i,
j) * t_diff_deviator(
M,
N,
k,
l));
516 t_flow_dir_dstrain(
i,
j,
k,
l) =
517 t_flow_dir_dstress(
i,
j,
m,
n) * t_D_Op(
m,
n,
k,
l);
523 commonDataPtr->resC.resize(nb_gauss_pts,
false);
524 commonDataPtr->resCdTau.resize(nb_gauss_pts,
false);
525 commonDataPtr->resCdStrain.resize(
size_symm, nb_gauss_pts,
false);
526 commonDataPtr->resCdPlasticStrain.resize(
size_symm, nb_gauss_pts,
false);
527 commonDataPtr->resFlow.resize(
size_symm, nb_gauss_pts,
false);
528 commonDataPtr->resFlowDtau.resize(
size_symm, nb_gauss_pts,
false);
534 commonDataPtr->resC.clear();
535 commonDataPtr->resCdTau.clear();
536 commonDataPtr->resCdStrain.clear();
537 commonDataPtr->resCdPlasticStrain.clear();
538 commonDataPtr->resFlow.clear();
539 commonDataPtr->resFlowDtau.clear();
540 commonDataPtr->resFlowDstrain.clear();
541 commonDataPtr->resFlowDstrainDot.clear();
543 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
544 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
545 auto t_res_c_dstrain =
546 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
547 auto t_res_c_plastic_strain =
548 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
549 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
550 auto t_res_flow_dtau =
551 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
552 auto t_res_flow_dstrain =
553 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
554 auto t_res_flow_dplastic_strain =
555 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
563 ++t_plastic_strain_dot;
568 ++t_res_c_plastic_strain;
571 ++t_res_flow_dstrain;
572 ++t_res_flow_dplastic_strain;
576 auto get_avtive_pts = [&]() {
577 int nb_points_avtive_on_elem = 0;
578 int nb_points_on_elem = 0;
580 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
581 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
582 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
583 auto t_plastic_strain_dot =
584 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
586 for (
auto &f : commonDataPtr->plasticSurface) {
589 eqiv, t_tau_dot, t_f,
597 ++nb_points_avtive_on_elem;
603 ++t_plastic_strain_dot;
613 nb_points += nb_points_on_elem;
614 if (nb_points_avtive_on_elem > 0) {
616 active_points += nb_points_avtive_on_elem;
617 if (nb_points_avtive_on_elem == nb_points_on_elem) {
622 if (nb_points_avtive_on_elem != nb_points_on_elem)
632 for (
auto &f : commonDataPtr->plasticSurface) {
636 t_diff_plastic_strain,
642 const auto d_sigma_y =
650 auto c =
constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
662 t_stress,
trace(t_stress),
669 t_flow_dir(
k,
l) = 1.5 * (t_dev_stress(
I,
J) * t_diff_deviator(
I,
J,
k,
l));
671 t_flow_dstrain(
i,
j) = t_flow(
k,
l) * t_D_Op(
k,
l,
i,
j);
673 auto get_res_c = [&]() {
return c; };
675 auto get_res_c_dstrain = [&](
auto &t_diff_res) {
676 t_diff_res(
i,
j) = c_f * t_flow_dstrain(
i,
j);
679 auto get_res_c_dplastic_strain = [&](
auto &t_diff_res) {
680 t_diff_res(
i,
j) = (this->getTSa() * c_equiv) * t_diff_eqiv(
i,
j);
681 t_diff_res(
k,
l) -= c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
684 auto get_res_c_dtau = [&]() {
685 return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
688 auto get_res_c_plastic_strain = [&](
auto &t_diff_res) {
689 t_diff_res(
k,
l) = -c_f * t_flow(
i,
j) * t_alpha_dir(
i,
j,
k,
l);
692 auto get_res_flow = [&](
auto &t_res_flow) {
693 const auto a = sigma_y;
694 const auto b = t_tau_dot;
695 t_res_flow(
k,
l) =
a * t_plastic_strain_dot(
k,
l) - b * t_flow_dir(
k,
l);
698 auto get_res_flow_dtau = [&](
auto &t_res_flow_dtau) {
699 const auto da = d_sigma_y;
700 const auto db = this->getTSa();
701 t_res_flow_dtau(
k,
l) =
702 da * t_plastic_strain_dot(
k,
l) - db * t_flow_dir(
k,
l);
705 auto get_res_flow_dstrain = [&](
auto &t_res_flow_dstrain) {
706 const auto b = t_tau_dot;
707 t_res_flow_dstrain(
m,
n,
k,
l) = -t_flow_dir_dstrain(
m,
n,
k,
l) * b;
710 auto get_res_flow_dplastic_strain = [&](
auto &t_res_flow_dplastic_strain) {
711 const auto a = sigma_y;
712 t_res_flow_dplastic_strain(
m,
n,
k,
l) =
713 (
a * this->getTSa()) * t_diff_plastic_strain(
m,
n,
k,
l);
714 const auto b = t_tau_dot;
715 t_res_flow_dplastic_strain(
m,
n,
i,
j) +=
716 (t_flow_dir_dstrain(
m,
n,
k,
l) * t_alpha_dir(
k,
l,
i,
j)) * b;
719 t_res_c = get_res_c();
720 get_res_flow(t_res_flow);
722 if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
723 t_res_c_dtau = get_res_c_dtau();
724 get_res_c_dstrain(t_res_c_dstrain);
725 get_res_c_dplastic_strain(t_res_c_plastic_strain);
726 get_res_flow_dtau(t_res_flow_dtau);
727 get_res_flow_dstrain(t_res_flow_dstrain);
728 get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
737template <
int DIM,
typename DomainEleOp>
740 boost::shared_ptr<CommonData> common_data_ptr,
741 boost::shared_ptr<MatrixDouble> mDPtr);
745 boost::shared_ptr<MatrixDouble>
mDPtr;
749template <
int DIM,
typename DomainEleOp>
751 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
752 boost::shared_ptr<MatrixDouble> m_D_ptr)
754 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
761template <
int DIM,
typename DomainEleOp>
772 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
773 commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
774 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
776 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
777 auto t_plastic_strain =
778 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
780 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
782 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
784 t_D(
i,
j,
k,
l) * (t_strain(
k,
l) - t_plastic_strain(
k,
l));
794template <
int DIM,
typename AssemblyDomainEleOp>
798 boost::shared_ptr<CommonData> common_data_ptr,
799 boost::shared_ptr<MatrixDouble> m_D_ptr);
800 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
804 boost::shared_ptr<MatrixDouble>
mDPtr;
807template <
int DIM,
typename AssemblyDomainEleOp>
810 boost::shared_ptr<CommonData> common_data_ptr,
811 boost::shared_ptr<MatrixDouble> m_D_ptr)
813 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
815template <
int DIM,
typename AssemblyDomainEleOp>
818 EntitiesFieldData::EntData &data) {
825 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
828 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
829 const auto nb_base_functions = data.getN().size2();
831 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
835 auto next = [&]() { ++t_res_flow; };
837 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
838 auto t_base = data.getFTensor0N();
840 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
841 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
845 t_rhs(
L) = alpha * (t_res_flow(
i,
j) * t_L(
i,
j,
L));
848 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
851 t_nf(
L) += t_base * t_rhs(
L);
855 for (; bb < nb_base_functions; ++bb)
862template <
typename AssemblyDomainEleOp>
866 boost::shared_ptr<CommonData> common_data_ptr,
867 boost::shared_ptr<MatrixDouble> m_D_ptr);
868 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
872 boost::shared_ptr<MatrixDouble>
mDPtr;
875template <
typename AssemblyDomainEleOp>
878 boost::shared_ptr<CommonData> common_data_ptr,
879 boost::shared_ptr<MatrixDouble> m_D_ptr)
881 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
883template <
typename AssemblyDomainEleOp>
886 EntitiesFieldData::EntData &data) {
889 const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
890 const size_t nb_base_functions = data.getN().size2();
892 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
894 auto next = [&]() { ++t_res_c; };
896 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
898 auto t_base = data.getFTensor0N();
899 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
900 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
902 const auto res = alpha * t_res_c;
907 nf[bb] += t_base * res;
910 for (; bb < nb_base_functions; ++bb)
917template <
int DIM,
typename AssemblyDomainEleOp>
921 const std::string row_field_name,
const std::string col_field_name,
922 boost::shared_ptr<CommonData> common_data_ptr,
923 boost::shared_ptr<MatrixDouble> m_D_ptr);
924 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
925 EntitiesFieldData::EntData &col_data);
929 boost::shared_ptr<MatrixDouble>
mDPtr;
932template <
int DIM,
typename AssemblyDomainEleOp>
935 const std::string row_field_name,
const std::string col_field_name,
936 boost::shared_ptr<CommonData> common_data_ptr,
937 boost::shared_ptr<MatrixDouble> m_D_ptr)
940 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
941 AssemblyDomainEleOp::sYmm =
false;
947 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
948 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
949 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
955 &mat(6 * rr + 0, 0), &mat(6 * rr + 0, 1), &mat(6 * rr + 0, 2),
956 &mat(6 * rr + 0, 3), &mat(6 * rr + 0, 4), &mat(6 * rr + 0, 5),
957 &mat(6 * rr + 1, 0), &mat(6 * rr + 1, 1), &mat(6 * rr + 1, 2),
958 &mat(6 * rr + 1, 3), &mat(6 * rr + 1, 4), &mat(6 * rr + 1, 5),
959 &mat(6 * rr + 2, 0), &mat(6 * rr + 2, 1), &mat(6 * rr + 2, 2),
960 &mat(6 * rr + 2, 3), &mat(6 * rr + 2, 4), &mat(6 * rr + 2, 5),
961 &mat(6 * rr + 3, 0), &mat(6 * rr + 3, 1), &mat(6 * rr + 3, 2),
962 &mat(6 * rr + 3, 3), &mat(6 * rr + 3, 4), &mat(6 * rr + 3, 5),
963 &mat(6 * rr + 4, 0), &mat(6 * rr + 4, 1), &mat(6 * rr + 4, 2),
964 &mat(6 * rr + 4, 3), &mat(6 * rr + 4, 4), &mat(6 * rr + 4, 5),
965 &mat(6 * rr + 5, 0), &mat(6 * rr + 5, 1), &mat(6 * rr + 5, 2),
966 &mat(6 * rr + 5, 3), &mat(6 * rr + 5, 4), &mat(6 * rr + 5, 5)};
969template <
int DIM,
typename AssemblyDomainEleOp>
972 EntitiesFieldData::EntData &row_data,
973 EntitiesFieldData::EntData &col_data) {
980 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
986 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
987 const auto nb_row_base_functions = row_data.getN().size2();
989 auto t_res_flow_dstrain =
990 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
991 auto t_res_flow_dplastic_strain =
992 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
996 ++t_res_flow_dstrain;
997 ++t_res_flow_dplastic_strain;
1000 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1001 auto t_row_base = row_data.getFTensor0N();
1002 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1003 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1008 alpha * (t_L(
i,
j, O) * ((t_res_flow_dplastic_strain(
i,
j,
k,
l) -
1009 t_res_flow_dstrain(
i,
j,
k,
l)) *
1017 auto t_col_base = col_data.getFTensor0N(gg, 0);
1019 t_mat(O,
L) += ((t_row_base * t_col_base) * t_res_mat(O,
L));
1027 for (; rr < nb_row_base_functions; ++rr)
1034template <
int DIM,
typename AssemblyDomainEleOp>
1038 const std::string row_field_name,
const std::string col_field_name,
1039 boost::shared_ptr<CommonData> common_data_ptr,
1040 boost::shared_ptr<MatrixDouble> m_D_ptr);
1041 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1042 EntitiesFieldData::EntData &col_data);
1049template <
int DIM,
typename AssemblyDomainEleOp>
1053 const std::string row_field_name,
const std::string col_field_name,
1054 boost::shared_ptr<CommonData> common_data_ptr,
1055 boost::shared_ptr<MatrixDouble> m_D_ptr);
1056 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1057 EntitiesFieldData::EntData &col_data);
1064template <
int DIM,
typename AssemblyDomainEleOp>
1067 const std::string row_field_name,
const std::string col_field_name,
1068 boost::shared_ptr<CommonData> common_data_ptr,
1069 boost::shared_ptr<MatrixDouble> m_D_ptr)
1072 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1073 AssemblyDomainEleOp::sYmm =
false;
1079 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
1085 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
1086 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
1089template <
int DIM,
typename AssemblyDomainEleOp>
1092 EntitiesFieldData::EntData &row_data,
1093 EntitiesFieldData::EntData &col_data) {
1098 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1101 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1102 const size_t nb_row_base_functions = row_data.getN().size2();
1105 const auto type = AssemblyDomainEleOp::getFEType();
1106 const auto nb_nodes = moab::CN::VerticesPerEntity(type);
1108 auto t_res_flow_dtau =
1109 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
1113 auto next = [&]() { ++t_res_flow_dtau; };
1115 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1116 auto t_row_base = row_data.getFTensor0N();
1117 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1118 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1121 t_res_vec(
L) = alpha * (t_res_flow_dtau(
i,
j) * t_L(
i,
j,
L));
1128 auto t_col_base = col_data.getFTensor0N(gg, 0);
1130 t_mat(
L) += t_row_base * t_col_base * t_res_vec(
L);
1136 for (; rr != nb_row_base_functions; ++rr)
1143template <
int DIM,
typename AssemblyDomainEleOp>
1147 const std::string row_field_name,
const std::string col_field_name,
1148 boost::shared_ptr<CommonData> common_data_ptr,
1149 boost::shared_ptr<MatrixDouble> mat_D_ptr);
1150 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1151 EntitiesFieldData::EntData &col_data);
1158template <
int DIM,
typename AssemblyDomainEleOp>
1161 const std::string row_field_name,
const std::string col_field_name,
1162 boost::shared_ptr<CommonData> common_data_ptr,
1163 boost::shared_ptr<MatrixDouble> m_D_ptr)
1166 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1167 AssemblyDomainEleOp::sYmm =
false;
1172 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
1177 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
1180template <
int DIM,
typename AssemblyDomainEleOp>
1183 EntitiesFieldData::EntData &row_data,
1184 EntitiesFieldData::EntData &col_data) {
1189 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
1192 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1193 const auto nb_row_base_functions = row_data.getN().size2();
1196 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdStrain);
1197 auto t_c_dplastic_strain =
1198 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->resCdPlasticStrain);
1202 ++t_c_dplastic_strain;
1207 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1208 auto t_row_base = row_data.getFTensor0N();
1209 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1210 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1215 t_L(
i,
j,
L) * (t_c_dplastic_strain(
i,
j) - t_c_dstrain(
i,
j));
1222 const auto row_base = alpha * t_row_base;
1223 auto t_col_base = col_data.getFTensor0N(gg, 0);
1225 t_mat(
L) += (row_base * t_col_base) * t_res_vec(
L);
1231 for (; rr != nb_row_base_functions; ++rr)
1238template <
typename AssemblyDomainEleOp>
1242 const std::string row_field_name,
const std::string col_field_name,
1243 boost::shared_ptr<CommonData> common_data_ptr);
1244 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1245 EntitiesFieldData::EntData &col_data);
1251template <
typename AssemblyDomainEleOp>
1254 const std::string row_field_name,
const std::string col_field_name,
1255 boost::shared_ptr<CommonData> common_data_ptr)
1258 commonDataPtr(common_data_ptr) {
1259 AssemblyDomainEleOp::sYmm =
false;
1262template <
typename AssemblyDomainEleOp>
1265 EntitiesFieldData::EntData &row_data,
1266 EntitiesFieldData::EntData &col_data) {
1269 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
1270 const auto nb_row_base_functions = row_data.getN().size2();
1272 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
1273 auto next = [&]() { ++t_res_c_dtau; };
1275 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
1276 auto t_row_base = row_data.getFTensor0N();
1277 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
1278 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
1281 const auto res = alpha * (t_res_c_dtau);
1287 auto t_col_base = col_data.getFTensor0N(gg, 0);
1289 *mat_ptr += t_row_base * t_col_base * res;
1295 for (; rr < nb_row_base_functions; ++rr)
Kronecker Delta class symmetric.
#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< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static auto get_mat_tensor_sym_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
auto diff_plastic_flow_dstress(long double f, FTensor::Tensor2_symmetric< T, DIM > &t_flow, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
static auto get_mat_tensor_sym_dscalar(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
double constrian_sign(double x)
double platsic_surface(FTensor::Tensor2_symmetric< double, 3 > &&t_stress_deviator)
auto get_mat_scalar_dtensor_sym(MatrixDouble &mat, FTensor::Number< 2 >)
auto symm_L_tensor(FTensor::Number< DIM >)
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
auto diff_constrain_dstress(double diff_constrain_df, FTensor::Tensor2_symmetric< T, DIM > &t_plastic_flow)
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]
double constrain_abs(double x)
auto diff_symmetrize(FTensor::Number< DIM >)
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 diff_plastic_flow_dstrain(FTensor::Ddg< T, DIM, DIM > &t_D, FTensor::Ddg< double, DIM, DIM > &&t_diff_plastic_flow_dstress)
auto plastic_flow(long double f, FTensor::Tensor2_symmetric< double, 3 > &&t_dev_stress, FTensor::Ddg< double, 3, DIM > &&t_diff_deviator)
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 diff_constrain_dstrain(T1 &t_D, T2 &&t_diff_constrain_dstress)
static auto get_mat_tensor_sym_dvector(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
[Lambda functions]
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
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 zeta
Viscous hardening.
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
const TSMethod::TSContext getTSCtx() const
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
static std::array< int, 5 > activityData
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr