Operator for linear form, usually to calculate values on right hand side.
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);
84 nb_integration_pts,
false);
85 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
87 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
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);
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.;
333 if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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)));
479 if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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)));
643 if (
ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
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");
Public Member Functions inherited from MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator