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 /
dt);
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 /
dt);
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 /
dt) * std::log1p(e);
276inline double w(
double eqiv,
double dot_tau,
double f,
double sigma_y,
278 return (1. /
cn1) * ((f - sigma_y) / sigma_Y) + dot_tau;
296inline double constraint(
double eqiv,
double dot_tau,
double f,
double sigma_y,
297 double abs_w,
double vis_H,
double sigma_Y) {
305 (
cn0 *
cn1 * ((dot_tau - eqiv)) +
307 cn1 * ((dot_tau) - (1. /
cn1) * (f - sigma_y) / sigma_Y - abs_w))
313 double vis_H,
double sigma_Y) {
314 return vis_H + (sigma_Y / 2) * (
cn0 *
cn1 +
cn1 * (1 - sign));
319 return (sigma_Y / 2) * (-
cn0 *
cn1);
326template <
typename T,
int DIM>
334 return t_diff_constrain_dstress;
337template <
typename T1,
typename T2,
int DIM>
344 t_diff_constrain_dstrain(
k,
l) =
345 t_diff_constrain_dstress(
i,
j) * t_D(
i,
j,
k,
l);
346 return t_diff_constrain_dstrain;
349template <
typename T,
int DIM>
354 constexpr double A = 2. / 3;
355 return std::sqrt(A * t_plastic_strain_dot(
i,
j) *
356 t_plastic_strain_dot(
i,
j)) +
357 std::numeric_limits<double>::epsilon();
360template <
typename T1,
typename T2,
typename T3,
int DIM>
362 T3 &t_diff_plastic_strain,
368 constexpr double A = 2. / 3;
370 t_diff_eqiv(
i,
j) = A * (t_plastic_strain_dot(
k,
l) / eqiv) *
371 t_diff_plastic_strain(
k,
l,
i,
j);
381 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 1, 0),
382 &mat(3 * rr + 1, 1), &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1)};
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)};
398template <
int DIM,
typename DomainEleOp>
402 boost::shared_ptr<CommonData> common_data_ptr);
403 MoFEMErrorCode doWork(
int side, EntityType type,
EntData &data);
409template <
int DIM,
typename DomainEleOp>
412 boost::shared_ptr<CommonData> common_data_ptr)
414 commonDataPtr(common_data_ptr) {
416 std::fill(&DomainEleOp::doEntities[MBEDGE],
417 &DomainEleOp::doEntities[MBMAXTYPE],
false);
420template <
int DIM,
typename DomainEleOp>
422 int side, EntityType type,
EntData &data) {
428 const size_t nb_gauss_pts = commonDataPtr->mStressPtr->size2();
430 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
431 auto t_plastic_strain =
432 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
434 commonDataPtr->plasticSurface.resize(nb_gauss_pts,
false);
435 commonDataPtr->plasticFlow.resize(
size_symm, nb_gauss_pts,
false);
436 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
438 auto ¶ms = commonDataPtr->blockParams;
440 for (
auto &f : commonDataPtr->plasticSurface) {
445 t_stress,
trace(t_stress),
458 t_flow(
i,
j) = t_flow_tmp(
i,
j);
468template <
int DIM,
typename DomainEleOp>
471 boost::shared_ptr<CommonData> common_data_ptr,
472 boost::shared_ptr<MatrixDouble> m_D_ptr);
473 MoFEMErrorCode doWork(
int side, EntityType type,
EntData &data);
477 boost::shared_ptr<MatrixDouble>
mDPtr;
480template <
int DIM,
typename DomainEleOp>
482 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
483 boost::shared_ptr<MatrixDouble> m_D_ptr)
485 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
487 std::fill(&DomainEleOp::doEntities[MBEDGE],
488 &DomainEleOp::doEntities[MBMAXTYPE],
false);
491template <
int DIM,
typename DomainEleOp>
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);
752template <
int DIM,
typename DomainEleOp>
759 boost::shared_ptr<CommonData> common_data_ptr,
760 boost::shared_ptr<MatrixDouble> mDPtr);
762 boost::shared_ptr<MatrixDouble> mDPtr);
767 boost::shared_ptr<MatrixDouble>
mDPtr;
771template <
int DIM,
typename DomainEleOp>
773 const std::string
field_name, boost::shared_ptr<CommonData> common_data_ptr,
774 boost::shared_ptr<MatrixDouble> m_D_ptr)
776 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
778 std::fill(&DomainEleOp::doEntities[MBEDGE],
779 &DomainEleOp::doEntities[MBMAXTYPE],
false);
782template <
int DIM,
typename DomainEleOp>
784 boost::shared_ptr<CommonData> common_data_ptr,
785 boost::shared_ptr<MatrixDouble> m_D_ptr)
787 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
790template <
int DIM,
typename DomainEleOp>
801 const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
802 commonDataPtr->mStressPtr->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts);
803 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
805 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStrainPtr));
806 auto t_plastic_strain =
807 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
809 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
811 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
813 t_D(
i,
j,
k,
l) * (t_strain(
k,
l) - t_plastic_strain(
k,
l));
823template <
int DIM,
typename AssemblyDomainEleOp>
827 boost::shared_ptr<CommonData> common_data_ptr,
828 boost::shared_ptr<MatrixDouble> m_D_ptr);
829 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
833 boost::shared_ptr<MatrixDouble>
mDPtr;
836template <
int DIM,
typename AssemblyDomainEleOp>
839 boost::shared_ptr<CommonData> common_data_ptr,
840 boost::shared_ptr<MatrixDouble> m_D_ptr)
842 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
844template <
int DIM,
typename AssemblyDomainEleOp>
847 EntitiesFieldData::EntData &data) {
854 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
857 const auto nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
858 const auto nb_base_functions = data.getN().size2();
860 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
864 auto next = [&]() { ++t_res_flow; };
866 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
867 auto t_base = data.getFTensor0N();
868 auto &nf = AssemblyDomainEleOp::locF;
869 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
870 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
874 t_rhs(L) = alpha * (t_res_flow(
i,
j) * t_L(
i,
j, L));
877 auto t_nf = getFTensor1FromArray<size_symm, size_symm>(nf);
879 for (; bb != AssemblyDomainEleOp::nbRows /
size_symm; ++bb) {
880 t_nf(L) += t_base * t_rhs(L);
884 for (; bb < nb_base_functions; ++bb)
891template <
typename AssemblyDomainEleOp>
895 boost::shared_ptr<CommonData> common_data_ptr,
896 boost::shared_ptr<MatrixDouble> m_D_ptr);
897 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
901 boost::shared_ptr<MatrixDouble>
mDPtr;
904template <
typename AssemblyDomainEleOp>
907 boost::shared_ptr<CommonData> common_data_ptr,
908 boost::shared_ptr<MatrixDouble> m_D_ptr)
910 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {}
912template <
typename AssemblyDomainEleOp>
915 EntitiesFieldData::EntData &data) {
918 const size_t nb_integration_pts = AssemblyDomainEleOp::getGaussPts().size2();
919 const size_t nb_base_functions = data.getN().size2();
921 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
923 auto next = [&]() { ++t_res_c; };
925 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
926 auto &nf = AssemblyDomainEleOp::locF;
927 auto t_base = data.getFTensor0N();
928 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
929 const double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
931 const auto res = alpha * t_res_c;
935 for (; bb != AssemblyDomainEleOp::nbRows; ++bb) {
936 nf[bb] += t_base * res;
939 for (; bb < nb_base_functions; ++bb)
946template <
int DIM,
typename AssemblyDomainEleOp>
950 const std::string row_field_name,
const std::string col_field_name,
951 boost::shared_ptr<CommonData> common_data_ptr,
952 boost::shared_ptr<MatrixDouble> m_D_ptr);
953 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
954 EntitiesFieldData::EntData &col_data);
958 boost::shared_ptr<MatrixDouble>
mDPtr;
961template <
int DIM,
typename AssemblyDomainEleOp>
964 const std::string row_field_name,
const std::string col_field_name,
965 boost::shared_ptr<CommonData> common_data_ptr,
966 boost::shared_ptr<MatrixDouble> m_D_ptr)
969 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
970 AssemblyDomainEleOp::sYmm =
false;
976 &mat(3 * rr + 0, 0), &mat(3 * rr + 0, 1), &mat(3 * rr + 0, 2),
977 &mat(3 * rr + 1, 0), &mat(3 * rr + 1, 1), &mat(3 * rr + 1, 2),
978 &mat(3 * rr + 2, 0), &mat(3 * rr + 2, 1), &mat(3 * rr + 2, 2)};
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)};
998template <
int DIM,
typename AssemblyDomainEleOp>
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)
1063template <
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);
1070 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1071 EntitiesFieldData::EntData &col_data);
1078template <
int DIM,
typename AssemblyDomainEleOp>
1082 const std::string row_field_name,
const std::string col_field_name,
1083 boost::shared_ptr<CommonData> common_data_ptr,
1084 boost::shared_ptr<MatrixDouble> m_D_ptr);
1085 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1086 EntitiesFieldData::EntData &col_data);
1093template <
int DIM,
typename AssemblyDomainEleOp>
1096 const std::string row_field_name,
const std::string col_field_name,
1097 boost::shared_ptr<CommonData> common_data_ptr,
1098 boost::shared_ptr<MatrixDouble> m_D_ptr)
1101 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1102 AssemblyDomainEleOp::sYmm =
false;
1108 &mat(3 * rr + 0, 0), &mat(3 * rr + 1, 0), &mat(3 * rr + 2, 0)};
1114 &mat(6 * rr + 0, 0), &mat(6 * rr + 1, 0), &mat(6 * rr + 2, 0),
1115 &mat(6 * rr + 3, 0), &mat(6 * rr + 4, 0), &mat(6 * rr + 5, 0)};
1118template <
int DIM,
typename AssemblyDomainEleOp>
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)
1169template <
int DIM,
typename AssemblyDomainEleOp>
1173 const std::string row_field_name,
const std::string col_field_name,
1174 boost::shared_ptr<CommonData> common_data_ptr,
1175 boost::shared_ptr<MatrixDouble> mat_D_ptr);
1176 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1177 EntitiesFieldData::EntData &col_data);
1184template <
int DIM,
typename AssemblyDomainEleOp>
1187 const std::string row_field_name,
const std::string col_field_name,
1188 boost::shared_ptr<CommonData> common_data_ptr,
1189 boost::shared_ptr<MatrixDouble> m_D_ptr)
1192 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
1193 AssemblyDomainEleOp::sYmm =
false;
1198 &mat(0, 0), &mat(0, 1), &mat(0, 2)};
1203 &mat(0, 0), &mat(0, 1), &mat(0, 2), &mat(0, 3), &mat(0, 4), &mat(0, 5)};
1206template <
int DIM,
typename AssemblyDomainEleOp>
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)
1264template <
typename AssemblyDomainEleOp>
1268 const std::string row_field_name,
const std::string col_field_name,
1269 boost::shared_ptr<CommonData> common_data_ptr);
1270 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1271 EntitiesFieldData::EntData &col_data);
1277template <
typename AssemblyDomainEleOp>
1280 const std::string row_field_name,
const std::string col_field_name,
1281 boost::shared_ptr<CommonData> common_data_ptr)
1284 commonDataPtr(common_data_ptr) {
1285 AssemblyDomainEleOp::sYmm =
false;
1288template <
typename AssemblyDomainEleOp>
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)
EntitiesFieldData::EntData EntData
Data on entities.
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< '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
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 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)
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]
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)
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 zeta
Viscous hardening.
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
auto kinematic_hardening_dplastic_strain(double C1_k)
constexpr auto field_name
FTensor::Index< 'm', 3 > m
Data on single entity (This is passed as argument to DataOperator::doWork)
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
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< MatrixDouble > mDPtr