13 #include <boost/math/constants/constants.hpp>
31 int nb_integration_pts = getGaussPts().size2();
33 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
34 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
39 auto t_eshelby_stress =
40 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
44 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
45 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
61 int nb_integration_pts = getGaussPts().size2();
77 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
78 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
79 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
81 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
82 dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts,
false);
83 dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 *
size_symm,
84 nb_integration_pts,
false);
85 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
87 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
88 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
89 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
90 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
92 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
93 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
94 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
95 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
96 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
97 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
99 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
100 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
103 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
104 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
106 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
107 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
108 auto t_levi_kirchhoff_domega =
109 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
110 auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
111 dataAtPts->leviKirchhoffdLogStreatchAtPts);
112 auto t_levi_kirchhoff_dP =
113 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
114 auto t_approx_P_adjoint_dstretch =
115 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
116 auto t_approx_P_adjoint_log_du =
117 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
118 auto t_approx_P_adjoint_log_du_dP =
119 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
120 auto t_approx_P_adjoint_log_du_domega =
121 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
122 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
123 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
125 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
126 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
127 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
128 auto t_log_stretch_total =
129 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
131 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
132 auto &nbUniq = dataAtPts->nbUniq;
137 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
138 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
139 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
141 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
149 ++t_levi_kirchhoff_domega;
150 ++t_levi_kirchhoff_dstreach;
151 ++t_levi_kirchhoff_dP;
152 ++t_approx_P_adjoint_dstretch;
153 ++t_approx_P_adjoint_log_du;
154 ++t_approx_P_adjoint_log_du_dP;
155 ++t_approx_P_adjoint_log_du_domega;
163 ++t_log_stretch_total;
174 auto bound_eig = [&](
auto &eig) {
176 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
177 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
178 for (
int ii = 0; ii != 3; ++ii) {
179 eig(ii) = std::min(std::max(zero, eig(ii)), large);
184 auto calculate_log_stretch = [&]() {
187 eigen_vec(
i,
j) = t_log_u(
i,
j);
189 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
193 t_nb_uniq = get_uniq_nb<3>(&eig(0));
197 t_eigen_vals(
i) = eig(
i);
198 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
201 auto get_t_diff_u = [&]() {
206 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
208 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
212 auto calculate_total_stretch = [&](
auto &t_h1) {
215 t_log_u2_h1(
i,
j) = 0;
216 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
221 t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
223 return std::make_pair(t_R_h1, t_inv_u_h1);
231 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
232 eigen_vec(
i,
j) = t_C_h1(
i,
j);
234 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
243 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
245 auto f_inv_sqrt = [](
auto e) {
return 1. / std::sqrt(e); };
249 t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
251 return std::make_pair(t_R_h1, t_inv_u_h1);
255 auto no_h1_loop = [&]() {
265 "rotSelector should be large");
268 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
276 auto t_Ldiff_u = calculate_log_stretch();
279 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
304 "rotationSelector not handled");
307 constexpr
double half_r = 1 / 2.;
308 constexpr
double half_l = 1 / 2.;
311 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
314 t_approx_P_adjoint_dstretch(
l,
k) =
315 (t_R(
i,
l) * t_approx_P(
i,
k) + t_R(
i,
k) * t_approx_P(
i,
l));
316 t_approx_P_adjoint_dstretch(
l,
k) /= 2.;
318 t_approx_P_adjoint_log_du(
L) =
319 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
322 t_levi_kirchhoff(
m) =
324 half_r * (t_diff_Rr(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
325 t_diff_Rr(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)))
329 half_l * (t_diff_Rl(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
330 t_diff_Rl(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)));
331 t_levi_kirchhoff(
m) /= 2.;
336 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
340 half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
343 "symmetrySelector not handled");
346 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
348 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
351 t_A(
k,
l,
m) = t_diff_Rr(
i,
l,
m) * t_approx_P(
i,
k) +
352 t_diff_Rr(
i,
k,
m) * t_approx_P(
i,
l);
355 t_B(
k,
l,
m) = t_diff_Rl(
i,
l,
m) * t_approx_P(
i,
k) +
356 t_diff_Rl(
i,
k,
m) * t_approx_P(
i,
l);
359 t_approx_P_adjoint_log_du_domega(
m,
L) =
360 half_r * (t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L)) +
361 half_l * (t_B(
k,
l,
m) * t_Ldiff_u(
k,
l,
L));
363 t_levi_kirchhoff_domega(
m,
n) =
365 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
366 t_diff_diff_Rr(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)))
371 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
372 t_diff_diff_Rl(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)));
373 t_levi_kirchhoff_domega(
m,
n) /= 2.;
382 auto large_loop = [&]() {
392 "rotSelector should be large");
395 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
409 "gradApproximator not handled");
413 auto t_Ldiff_u = calculate_log_stretch();
415 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
418 t_h_u(
l,
k) = t_u(
l, o) * t_h1(o,
k);
420 t_Ldiff_h_u(
l,
k,
L) = t_Ldiff_u(
l, o,
L) * t_h1(o,
k);
447 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
450 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
455 "rotationSelector not handled");
458 constexpr
double half_r = 1 / 2.;
459 constexpr
double half_l = 1 / 2.;
462 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
465 t_approx_P_adjoint_dstretch(
l, o) =
466 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
467 t_approx_P_adjoint_log_du(
L) =
468 t_approx_P_adjoint_dstretch(
l, o) * t_Ldiff_u(
l, o,
L);
471 t_levi_kirchhoff(
m) =
473 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)))
477 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)));
482 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
486 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
488 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
491 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
493 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
494 t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
498 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
500 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
502 t_approx_P_adjoint_log_du_domega(
m,
L) =
503 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
507 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
510 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
511 t_approx_P_adjoint_log_du_domega(
m,
L) =
512 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
515 t_levi_kirchhoff_dstreach(
m,
L) =
517 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
521 half_l * (t_diff_Rl(
i,
l,
m) *
522 (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
524 t_levi_kirchhoff_dP(
m,
i,
k) =
526 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
530 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
532 t_levi_kirchhoff_domega(
m,
n) =
534 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
539 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
548 auto moderate_loop = [&]() {
558 "rotSelector should be large");
561 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
572 "gradApproximator not handled");
576 auto t_Ldiff_u = calculate_log_stretch();
578 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
581 t_h_u(
l,
k) = (t_symm_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
583 t_L_h(
l,
k,
L) = t_L(
l, o,
L) * t_h1(o,
k);
610 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
613 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
618 "rotationSelector not handled");
621 constexpr
double half_r = 1 / 2.;
622 constexpr
double half_l = 1 / 2.;
625 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
628 t_approx_P_adjoint_dstretch(
l, o) =
629 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
631 t_approx_P_adjoint_log_du(
L) =
632 t_approx_P_adjoint_dstretch(
l, o) * t_L(
l, o,
L);
635 t_levi_kirchhoff(
m) =
637 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)))
641 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k))*t_approx_P(
i,
k)));
646 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
650 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
652 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
655 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
657 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
661 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_L_h(
l,
k,
L);
663 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_L_h(
l,
k,
L);
664 t_approx_P_adjoint_log_du_domega(
m,
L) =
665 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
669 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
672 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_L_h(
l,
k,
L);
673 t_approx_P_adjoint_log_du_domega(
m,
L) =
674 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
677 t_levi_kirchhoff_dstreach(
m,
L) =
678 half_r * (t_diff_Rr(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)))
682 half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
684 t_levi_kirchhoff_dP(
m,
i,
k) =
686 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
690 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
692 t_levi_kirchhoff_domega(
m,
n) =
694 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
699 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
708 auto small_loop = [&]() {
715 "rotSelector should be small");
718 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
727 "gradApproximator not handled");
734 t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
736 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
737 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
739 t_log_u2_h1(
i,
j) = 0;
740 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
747 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
750 t_approx_P_adjoint_dstretch(
i,
j) = t_approx_P(
i,
j);
751 t_approx_P_adjoint_log_du(
L) =
752 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
753 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
754 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
758 t_levi_kirchhoff_dstreach(
m,
L) = 0;
760 t_levi_kirchhoff_domega(
m,
n) = 0;
786 "gradApproximator not handled");
799 auto n_total_side_eles = getLoopSize();
800 auto N_InLoop = getNinTheLoop();
801 auto sensee = getSkeletonSense();
802 auto nb_gauss_pts = getGaussPts().size2();
803 auto t_normal = getFTensor1NormalsAtGaussPts();
806 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
808 dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
809 dataAtPts->tractionAtPts.clear();
811 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
813 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
814 t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
826 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
832 int nb_integration_pts = getGaussPts().size2();
833 auto t_w = getFTensor0IntegrationWeight();
834 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
835 auto t_coords = getFTensor1CoordsAtGaussPts();
836 auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
844 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
845 loc_reaction_forces(
i) += (t_traction(
i))*t_w * getMeasure();
846 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
848 loc_moment_forces(
i) +=
849 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
850 t_traction(
k) * t_w * getMeasure();
857 reactionVec[0] += loc_reaction_forces(0);
858 reactionVec[1] += loc_reaction_forces(1);
859 reactionVec[2] += loc_reaction_forces(2);
860 reactionVec[3] += loc_moment_forces(0);
861 reactionVec[4] += loc_moment_forces(1);
862 reactionVec[5] += loc_moment_forces(2);
870 int nb_integration_pts = data.
getN().size1();
871 auto v = getVolume();
872 auto t_w = getFTensor0IntegrationWeight();
873 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
874 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
875 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
876 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
877 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
878 dataAtPts->wL2DotDotAtPts.clear();
880 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
882 auto piola_scale = dataAtPts->piolaScale;
883 auto alpha_w = alphaW / piola_scale;
884 auto alpha_rho = alphaRho / piola_scale;
886 int nb_base_functions = data.
getN().size2();
890 auto get_ftensor1 = [](
auto &
v) {
895 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
897 auto t_nf = get_ftensor1(nF);
899 for (; bb != nb_dofs / 3; ++bb) {
900 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
901 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
902 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
906 for (; bb != nb_base_functions; ++bb)
920 int nb_integration_pts = getGaussPts().size2();
921 auto v = getVolume();
922 auto t_w = getFTensor0IntegrationWeight();
923 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
924 auto t_omega_grad_dot =
925 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
926 int nb_base_functions = data.
getN().size2();
932 auto get_ftensor1 = [](
auto &
v) {
936 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
938 auto t_nf = get_ftensor1(nF);
940 for (; bb != nb_dofs / 3; ++bb) {
941 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
943 (
a * alphaOmega) * (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
948 for (; bb != nb_base_functions; ++bb) {
962 int nb_integration_pts = data.
getN().size1();
963 auto v = getVolume();
964 auto t_w = getFTensor0IntegrationWeight();
966 int nb_base_functions = data.
getN().size2() / 3;
974 auto get_ftensor1 = [](
auto &
v) {
981 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
982 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
984 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
986 auto t_nf = get_ftensor1(nF);
991 for (; bb != nb_dofs / 3; ++bb) {
992 t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
993 t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
994 t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
999 for (; bb != nb_base_functions; ++bb)
1009 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1011 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1013 auto t_nf = get_ftensor1(nF);
1021 for (; bb != nb_dofs / 3; ++bb) {
1022 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1027 for (; bb != nb_base_functions; ++bb)
1041 int nb_integration_pts = data.
getN().size1();
1042 auto v = getVolume();
1043 auto t_w = getFTensor0IntegrationWeight();
1045 int nb_base_functions = data.
getN().size2() / 9;
1053 auto get_ftensor0 = [](
auto &
v) {
1059 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1060 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1062 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1064 auto t_nf = get_ftensor0(nF);
1069 for (; bb != nb_dofs; ++bb) {
1070 t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
1071 t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
1075 for (; bb != nb_base_functions; ++bb) {
1085 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1087 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1089 auto t_nf = get_ftensor0(nF);
1093 t_residuum(
i,
j) = t_h(
i,
j);
1096 for (; bb != nb_dofs; ++bb) {
1097 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1101 for (; bb != nb_base_functions; ++bb) {
1115 int nb_integration_pts = data.
getN().size1();
1116 auto v = getVolume();
1117 auto t_w = getFTensor0IntegrationWeight();
1118 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1119 int nb_base_functions = data.
getN().size2() / 3;
1122 auto get_ftensor1 = [](
auto &
v) {
1127 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1129 auto t_nf = get_ftensor1(nF);
1131 for (; bb != nb_dofs / 3; ++bb) {
1132 double div_row_base = t_row_diff_base_fun(
i,
i);
1133 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1135 ++t_row_diff_base_fun;
1137 for (; bb != nb_base_functions; ++bb) {
1138 ++t_row_diff_base_fun;
1147 template <AssemblyType A>
1153 for (
auto &bc : (*bcDispPtr)) {
1155 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1157 int nb_integration_pts = OP::getGaussPts().size2();
1158 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1159 auto t_w = OP::getFTensor0IntegrationWeight();
1160 int nb_base_functions = data.
getN().size2() / 3;
1167 for (
auto &sm : scalingMethodsVec) {
1171 scale *= sm->getScale(OP::getFEMethod()->ts_t);
1179 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1180 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1182 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1184 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1188 for (; bb != nb_base_functions; ++bb)
1200 return OP::iNtegrate(data);
1203 template <AssemblyType A>
1214 for (
auto &bc : (*bcRotPtr)) {
1216 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1218 int nb_integration_pts = OP::getGaussPts().size2();
1219 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1220 auto t_w = OP::getFTensor0IntegrationWeight();
1222 int nb_base_functions = data.
getN().size2() / 3;
1225 auto get_ftensor1 = [](
auto &
v) {
1233 double theta = bc.theta;
1237 theta *= OP::getFEMethod()->ts_t;
1241 const double a = sqrt(t_normal(
i) * t_normal(
i));
1242 t_omega(
i) = t_normal(
i) * (theta /
a);
1247 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1249 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1251 t_delta(
i) = t_center(
i) - t_coords(
i);
1253 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1255 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1257 for (; bb != nb_dofs / 3; ++bb) {
1258 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1262 for (; bb != nb_base_functions; ++bb)
1275 return OP::iNtegrate(data);
1284 int nb_integration_pts = getGaussPts().size2();
1285 int nb_base_functions = data.
getN().size2();
1287 double ts_t = getFEMethod()->ts_t;
1292 double time_scale = 1;
1293 for (
auto &sm : scalingMethodsVec) {
1294 time_scale *= sm->getScale(ts_t);
1298 if (this->locF.size() != nb_dofs)
1300 "Size of locF %d != nb_dofs %d", this->locF.size(), nb_dofs);
1303 auto integrate_rhs = [&](
auto &bc,
auto calc_tau) {
1306 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1308 auto t_w = getFTensor0IntegrationWeight();
1309 auto t_coords = getFTensor1CoordsAtGaussPts();
1311 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1313 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1315 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1316 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1318 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1319 t_f(
i) -= time_scale * t_w * t_row_base * tau * (t_val(
i) *
scale);
1324 for (; rr != nb_base_functions; ++rr)
1331 this->locF *= 2. * getMeasure();
1337 for (
auto &bc : *(bcData)) {
1338 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1343 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1347 return -y * (y - 1) / 0.25;
1349 CHKERR integrate_rhs(bc, calc_tau);
1351 CHKERR integrate_rhs(bc, [](
double,
double,
double) {
return 1; });
1362 int nb_integration_pts = row_data.
getN().size1();
1363 int row_nb_dofs = row_data.
getIndices().size();
1364 int col_nb_dofs = col_data.
getIndices().size();
1367 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2));
1370 auto v = getVolume();
1371 auto t_w = getFTensor0IntegrationWeight();
1372 int row_nb_base_functions = row_data.
getN().size2();
1374 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1377 for (; rr != row_nb_dofs / 3; ++rr) {
1379 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1380 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1381 double div_col_base = t_col_diff_base_fun(
i,
i);
1382 t_m(
i) -=
a * t_row_base_fun * div_col_base;
1384 ++t_col_diff_base_fun;
1388 for (; rr != row_nb_base_functions; ++rr)
1399 if (alphaW < std::numeric_limits<double>::epsilon() &&
1400 alphaRho < std::numeric_limits<double>::epsilon())
1403 const int nb_integration_pts = row_data.
getN().size1();
1404 const int row_nb_dofs = row_data.
getIndices().size();
1407 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2)
1413 auto v = getVolume();
1414 auto t_w = getFTensor0IntegrationWeight();
1416 auto piola_scale = dataAtPts->piolaScale;
1417 auto alpha_w = alphaW / piola_scale;
1418 auto alpha_rho = alphaRho / piola_scale;
1420 int row_nb_base_functions = row_data.
getN().size2();
1423 double ts_scale = alpha_w * getTSa();
1424 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1425 ts_scale += alpha_rho * getTSaa();
1427 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1428 double a =
v * t_w * ts_scale;
1431 for (; rr != row_nb_dofs / 3; ++rr) {
1434 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1435 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1436 const double b =
a * t_row_base_fun * t_col_base_fun;
1445 for (; rr != row_nb_base_functions; ++rr)
1464 int nb_integration_pts = row_data.
getN().size1();
1465 int row_nb_dofs = row_data.
getIndices().size();
1466 int col_nb_dofs = col_data.
getIndices().size();
1470 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1472 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1474 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
1476 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
1478 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
1480 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2));
1483 auto v = getVolume();
1484 auto t_w = getFTensor0IntegrationWeight();
1486 int row_nb_base_functions = row_data.
getN().size2();
1489 auto t_approx_P_adjoint_log_du_dP =
1490 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1492 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1495 for (; rr != row_nb_dofs / 6; ++rr) {
1498 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1500 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1502 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
1510 for (; rr != row_nb_base_functions; ++rr)
1513 ++t_approx_P_adjoint_log_du_dP;
1529 int nb_integration_pts = row_data.
getN().size1();
1530 int row_nb_dofs = row_data.
getIndices().size();
1531 int col_nb_dofs = col_data.
getIndices().size();
1534 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c), &
m(
r + 3,
c), &
m(
r + 4,
c),
1538 auto v = getVolume();
1539 auto t_w = getFTensor0IntegrationWeight();
1542 int row_nb_base_functions = row_data.
getN().size2();
1544 auto t_approx_P_adjoint_log_du_dP =
1545 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1547 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1550 for (; rr != row_nb_dofs / 6; ++rr) {
1551 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1552 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1553 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1555 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
1562 for (; rr != row_nb_base_functions; ++rr)
1565 ++t_approx_P_adjoint_log_du_dP;
1577 int row_nb_dofs = row_data.
getIndices().size();
1578 int col_nb_dofs = col_data.
getIndices().size();
1582 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1584 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1586 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
1588 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
1590 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
1592 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2)
1602 auto v = getVolume();
1603 auto t_w = getFTensor0IntegrationWeight();
1604 auto t_approx_P_adjoint_log_du_domega =
1605 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
1607 int row_nb_base_functions = row_data.
getN().size2();
1610 int nb_integration_pts = row_data.
getN().size1();
1612 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1616 for (; rr != row_nb_dofs / 6; ++rr) {
1618 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1619 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1620 double v =
a * t_row_base_fun * t_col_base_fun;
1621 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
1628 for (; rr != row_nb_base_functions; ++rr)
1632 ++t_approx_P_adjoint_log_du_domega;
1641 int nb_integration_pts = getGaussPts().size2();
1642 int row_nb_dofs = row_data.
getIndices().size();
1643 int col_nb_dofs = col_data.
getIndices().size();
1648 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
1649 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
1651 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
1652 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
1654 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
1655 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5)
1663 auto v = getVolume();
1664 auto t_w = getFTensor0IntegrationWeight();
1665 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
1666 dataAtPts->leviKirchhoffdLogStreatchAtPts);
1667 int row_nb_base_functions = row_data.
getN().size2();
1669 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1672 for (; rr != row_nb_dofs / 3; ++rr) {
1673 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1674 const double b =
a * t_row_base_fun;
1676 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
1677 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
1683 for (; rr != row_nb_base_functions; ++rr) {
1687 ++t_levi_kirchhoff_du;
1703 int nb_integration_pts = getGaussPts().size2();
1704 int row_nb_dofs = row_data.
getIndices().size();
1705 int col_nb_dofs = col_data.
getIndices().size();
1709 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1711 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1713 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1718 auto v = getVolume();
1719 auto t_w = getFTensor0IntegrationWeight();
1721 int row_nb_base_functions = row_data.
getN().size2();
1723 auto t_levi_kirchhoff_dP =
1724 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1726 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1729 for (; rr != row_nb_dofs / 3; ++rr) {
1730 double b =
a * t_row_base_fun;
1732 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1733 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1734 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
1740 for (; rr != row_nb_base_functions; ++rr) {
1745 ++t_levi_kirchhoff_dP;
1753 int nb_integration_pts = getGaussPts().size2();
1754 int row_nb_dofs = row_data.
getIndices().size();
1755 int col_nb_dofs = col_data.
getIndices().size();
1759 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c));
1766 auto v = getVolume();
1767 auto t_w = getFTensor0IntegrationWeight();
1768 auto t_levi_kirchoff_dP =
1769 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1771 int row_nb_base_functions = row_data.
getN().size2();
1774 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1777 for (; rr != row_nb_dofs / 3; ++rr) {
1778 double b =
a * t_row_base_fun;
1779 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1780 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1781 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1782 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
1789 for (; rr != row_nb_base_functions; ++rr) {
1793 ++t_levi_kirchoff_dP;
1801 int nb_integration_pts = getGaussPts().size2();
1802 int row_nb_dofs = row_data.
getIndices().size();
1803 int col_nb_dofs = col_data.
getIndices().size();
1807 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1809 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1811 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1824 auto v = getVolume();
1825 auto ts_a = getTSa();
1826 auto t_w = getFTensor0IntegrationWeight();
1827 auto t_levi_kirchhoff_domega =
1828 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
1829 int row_nb_base_functions = row_data.
getN().size2();
1832 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1834 double c =
a * alphaOmega * ts_a;
1836 for (; rr != row_nb_dofs / 3; ++rr) {
1837 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1838 const double b =
a * t_row_base_fun;
1841 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1842 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
1843 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
1851 for (; rr != row_nb_base_functions; ++rr) {
1856 ++t_levi_kirchhoff_domega;
1864 if (dataAtPts->matInvD.size1() ==
size_symm &&
1865 dataAtPts->matInvD.size2() ==
size_symm) {
1866 return integrateImpl<0>(row_data, col_data);
1868 return integrateImpl<size_symm>(row_data, col_data);
1881 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1883 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1885 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1890 int nb_integration_pts = getGaussPts().size2();
1891 int row_nb_dofs = row_data.
getIndices().size();
1892 int col_nb_dofs = col_data.
getIndices().size();
1894 auto v = getVolume();
1895 auto t_w = getFTensor0IntegrationWeight();
1896 int row_nb_base_functions = row_data.
getN().size2() / 3;
1903 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1904 &*dataAtPts->matInvD.data().begin());
1907 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1911 for (; rr != row_nb_dofs / 3; ++rr) {
1913 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1914 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1915 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
1923 for (; rr != row_nb_base_functions; ++rr)
1936 if (dataAtPts->matInvD.size1() ==
size_symm &&
1937 dataAtPts->matInvD.size2() ==
size_symm) {
1938 return integrateImpl<0>(row_data, col_data);
1940 return integrateImpl<size_symm>(row_data, col_data);
1951 int nb_integration_pts = getGaussPts().size2();
1952 int row_nb_dofs = row_data.
getIndices().size();
1953 int col_nb_dofs = col_data.
getIndices().size();
1955 auto v = getVolume();
1956 auto t_w = getFTensor0IntegrationWeight();
1957 int row_nb_base_functions = row_data.
getN().size2() / 9;
1964 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1965 &*dataAtPts->matInvD.data().begin());
1968 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1972 for (; rr != row_nb_dofs; ++rr) {
1974 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1976 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
1983 for (; rr != row_nb_base_functions; ++rr)
1994 if (dataAtPts->matInvD.size1() ==
size_symm &&
1995 dataAtPts->matInvD.size2() ==
size_symm) {
1996 return integrateImpl<0>(row_data, col_data);
1998 return integrateImpl<size_symm>(row_data, col_data);
2012 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2)
2017 int nb_integration_pts = getGaussPts().size2();
2018 int row_nb_dofs = row_data.
getIndices().size();
2019 int col_nb_dofs = col_data.
getIndices().size();
2021 auto v = getVolume();
2022 auto t_w = getFTensor0IntegrationWeight();
2023 int row_nb_base_functions = row_data.
getN().size2() / 9;
2032 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2033 &*dataAtPts->matInvD.data().begin());
2036 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2039 auto t_m = get_ftensor1(
K, 0, 0);
2042 for (; rr != row_nb_dofs; ++rr) {
2044 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2045 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2053 for (; rr != row_nb_base_functions; ++rr)
2072 int nb_integration_pts = row_data.
getN().size1();
2073 int row_nb_dofs = row_data.
getIndices().size();
2074 int col_nb_dofs = col_data.
getIndices().size();
2079 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2081 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2083 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
2088 auto v = getVolume();
2089 auto t_w = getFTensor0IntegrationWeight();
2090 int row_nb_base_functions = row_data.
getN().size2() / 3;
2093 const double ts_a = getTSa();
2095 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2096 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2100 for (; rr != row_nb_dofs / 3; ++rr) {
2103 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2106 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2107 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2108 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
2116 for (; rr != row_nb_base_functions; ++rr)
2136 int nb_integration_pts = row_data.
getN().size1();
2137 int row_nb_dofs = row_data.
getIndices().size();
2138 int col_nb_dofs = col_data.
getIndices().size();
2142 &
m(
r,
c + 0), &
m(
r,
c + 1), &
m(
r,
c + 2));
2145 auto v = getVolume();
2146 auto t_w = getFTensor0IntegrationWeight();
2147 int row_nb_base_functions = row_data.
getN().size2() / 9;
2150 const double ts_a = getTSa();
2152 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2153 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2157 for (; rr != row_nb_dofs; ++rr) {
2160 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
2163 auto t_m = get_ftensor2(
K, rr, 0);
2164 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2165 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
2173 for (; rr != row_nb_base_functions; ++rr)
2186 if (tagSense != getSkeletonSense())
2189 auto get_tag = [&](
auto name) {
2190 auto &mob = getPtrFE()->mField.get_moab();
2196 auto get_tag_value = [&](
auto &&tag,
int dim) {
2197 auto &mob = getPtrFE()->mField.get_moab();
2198 auto face = getSidePtrFE()->getFEEntityHandle();
2199 std::vector<double> value(dim);
2200 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
2204 auto create_tag = [
this](
const std::string tag_name,
const int size) {
2205 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2207 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2208 th, MB_TAG_CREAT | MB_TAG_SPARSE,
2213 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
2214 Tag th_detF = create_tag(
"detF", 1);
2215 Tag th_traction = create_tag(
"traction", 3);
2216 Tag th_disp_error = create_tag(
"DisplacementError", 1);
2218 Tag th_energy = create_tag(
"Energy", 1);
2220 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
2221 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2222 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2224 auto t_normal = getFTensor1NormalsAtGaussPts();
2225 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
2235 if (dataAtPts->energyAtPts.size() == 0) {
2237 dataAtPts->energyAtPts.resize(getGaussPts().size2());
2238 dataAtPts->energyAtPts.clear();
2247 auto set_float_precision = [](
const double x) {
2248 if (std::abs(x) < std::numeric_limits<float>::epsilon())
2255 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
2257 v = set_float_precision(
v);
2258 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
2265 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2268 for (
auto &
a :
v.data())
2269 a = set_float_precision(
a);
2270 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
2271 &*
v.data().begin());
2279 &
m(0, 0), &
m(0, 1), &
m(0, 2),
2281 &
m(1, 0), &
m(1, 1), &
m(1, 2),
2283 &
m(2, 0), &
m(2, 1), &
m(2, 2));
2285 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2287 t_m(
i,
j) = t_d(
i,
j);
2288 for (
auto &
v :
m.data())
2289 v = set_float_precision(
v);
2290 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
2291 &*
m.data().begin());
2295 const auto nb_gauss_pts = getGaussPts().size2();
2296 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
2299 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
2301 CHKERR save_vec_tag(th_traction, t_traction, gg);
2303 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
2304 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2305 CHKERR save_scal_tag(th_energy, t_energy, gg);
2309 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
2310 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
2311 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2320 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
2321 std::vector<FieldSpace> spaces, std::string geom_field_name,
2322 boost::shared_ptr<Range> crack_front_edges_ptr) {
2325 constexpr
bool scale_l2 =
false;
2335 boost::shared_ptr<Range> edges_ptr)
2344 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
2347 return OP::doWork(side,
type, data);
2352 boost::shared_ptr<Range> edgesPtr;
2355 auto jac = boost::make_shared<MatrixDouble>();
2356 auto det = boost::make_shared<VectorDouble>();
2359 ? crack_front_edges_ptr
2360 : boost::make_shared<Range>()));
2373 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
2374 std::vector<FieldSpace> spaces, std::string geom_field_name,
2375 boost::shared_ptr<Range> crack_front_edges_ptr) {
2379 auto jac = boost::make_shared<MatrixDouble>();
2380 auto det = boost::make_shared<VectorDouble>();
2382 geom_field_name, jac));
2388 constexpr
bool scale_l2_ainsworth_legendre_base =
false;
2390 if (scale_l2_ainsworth_legendre_base) {
2398 boost::shared_ptr<MatrixDouble> jac,
2399 boost::shared_ptr<Range> edges_ptr)
2408 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
2411 return OP::doWork(side,
type, data);
2416 boost::shared_ptr<Range> edgesPtr;
2419 auto jac = boost::make_shared<MatrixDouble>();
2420 auto det = boost::make_shared<VectorDouble>();
2422 geom_field_name, jac,
2424 : boost::make_shared<Range>()));
2431 nullptr,
nullptr,
nullptr);
2450 const auto nb_integration_pts = getGaussPts().size2();
2452 auto t_w = getFTensor0IntegrationWeight();
2453 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2454 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2456 auto t_var_omega = getFTensor1FromMat<3>(dataAtPts->varRotAxis);
2457 auto t_var_log_u = getFTensor2SymmetricFromMat<3>(dataAtPts->varLogStreach);
2458 auto t_var_P = getFTensor2FromMat<3, 3>(dataAtPts->varPiola);
2460 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2462 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
2464 auto t_approx_P_adjoint_log_du =
2465 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
2467 double var_element_energy = 0.;
2469 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
2472 t_var_L_u(
L) = t_L(
J,
K,
L) * t_var_log_u(
J,
K);
2473 auto var_energy = t_P(
i,
J) * (t_h_domega(
i,
J,
j) * t_var_omega(
j) +
2474 t_h_dlog_u(
i,
J,
L) * t_var_L_u(
L));
2475 var_element_energy += t_w * var_energy;
2476 auto var_complementary = t_var_P(
i,
J) * t_h(
i,
J);
2477 var_element_energy += t_w * var_complementary;
2490 ++t_approx_P_adjoint_log_du;
2493 var_element_energy *= getMeasure();
2495 auto get_tag = [&]() {
2496 auto &mob = getPtrFE()->mField.get_moab();
2498 double def_val[] = {0.};
2499 CHK_MOAB_THROW(mob.tag_get_handle(
"ReleaseEnergy", 1, MB_TYPE_DOUBLE, tag,
2500 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
2505 auto set_tag = [&](
auto &&tag,
auto &energy) {
2506 auto &mob = getPtrFE()->mField.get_moab();
2507 auto tet = getPtrFE()->getFEEntityHandle();
2508 CHK_MOAB_THROW(mob.tag_set_data(tag, &tet, 1, &energy),
"set tag");
2511 set_tag(get_tag(), var_element_energy);
2513 *releaseEnergyPtr += var_element_energy;
