13#include <boost/math/constants/constants.hpp>
22#ifdef ENABLE_PYTHON_BINDING
23#include <boost/python.hpp>
24#include <boost/python/def.hpp>
25#include <boost/python/numpy.hpp>
26namespace bp = boost::python;
27namespace np = boost::python::numpy;
39 int nb_integration_pts = getGaussPts().size2();
41 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
42 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
47 auto t_eshelby_stress =
48 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
52 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
53 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
69 int nb_integration_pts = getGaussPts().size2();
85 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
86 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
87 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
89 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
90 dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts,
false);
91 dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 *
size_symm,
92 nb_integration_pts,
false);
93 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
95 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
96 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
97 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
98 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
100 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
101 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
102 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
103 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
104 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
105 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
106 dataAtPts->eigenValsC.resize(3, nb_integration_pts,
false);
107 dataAtPts->eigenVecsC.resize(9, nb_integration_pts,
false);
108 dataAtPts->nbUniqC.resize(nb_integration_pts,
false);
110 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
111 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
113 dataAtPts->internalStressAtPts.resize(9, nb_integration_pts,
false);
114 dataAtPts->internalStressAtPts.clear();
117 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
118 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
120 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
121 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
122 auto t_levi_kirchhoff_domega =
123 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
124 auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
125 dataAtPts->leviKirchhoffdLogStreatchAtPts);
126 auto t_levi_kirchhoff_dP =
127 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
128 auto t_approx_P_adjoint__dstretch =
129 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
130 auto t_approx_P_adjoint_log_du =
131 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
132 auto t_approx_P_adjoint_log_du_dP =
133 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
134 auto t_approx_P_adjoint_log_du_domega =
135 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
136 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
137 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
139 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
140 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
141 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
142 auto &nbUniq = dataAtPts->nbUniq;
145 auto t_eigen_vals_C = getFTensor1FromMat<3>(dataAtPts->eigenValsC);
146 auto t_eigen_vecs_C = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecsC);
147 auto &nbUniqC = dataAtPts->nbUniqC;
151 auto t_log_stretch_total =
152 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
154 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
157 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
158 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
159 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
161 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
169 ++t_levi_kirchhoff_domega;
170 ++t_levi_kirchhoff_dstreach;
171 ++t_levi_kirchhoff_dP;
172 ++t_approx_P_adjoint__dstretch;
173 ++t_approx_P_adjoint_log_du;
174 ++t_approx_P_adjoint_log_du_dP;
175 ++t_approx_P_adjoint_log_du_domega;
186 ++t_log_stretch_total;
197 auto bound_eig = [&](
auto &eig) {
199 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
200 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
201 for (
int ii = 0; ii != 3; ++ii) {
202 eig(ii) = std::min(std::max(zero, eig(ii)), large);
207 auto calculate_log_stretch = [&]() {
210 eigen_vec(
i,
j) = t_log_u(
i,
j);
212 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
216 t_nb_uniq = get_uniq_nb<3>(&eig(0));
220 t_eigen_vals(
i) = eig(
i);
221 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
224 auto get_t_diff_u = [&]() {
229 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
231 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
235 auto calculate_total_stretch = [&](
auto &t_h1) {
238 t_log_u2_h1(
i,
j) = 0;
239 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
244 t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
246 return std::make_pair(t_R_h1, t_inv_u_h1);
254 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
255 eigen_vec(
i,
j) = t_C_h1(
i,
j);
257 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
260 t_nb_uniq_C = get_uniq_nb<3>(&eig(0));
261 if (t_nb_uniq_C < 3) {
267 t_eigen_vals_C(
i) = eig(
i);
268 t_eigen_vecs_C(
i,
j) = eigen_vec(
i,
j);
272 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
274 auto f_inv_sqrt = [](
auto e) {
return 1. / std::sqrt(e); };
278 t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
280 return std::make_pair(t_R_h1, t_inv_u_h1);
284 auto no_h1_loop = [&]() {
294 "rotSelector should be large");
297 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
305 auto t_Ldiff_u = calculate_log_stretch();
308 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
324 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
325 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
327 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
328 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
333 "rotationSelector not handled");
336 constexpr double half_r = 1 / 2.;
337 constexpr double half_l = 1 / 2.;
340 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
343 t_approx_P_adjoint__dstretch(
l,
k) =
344 (t_R(
i,
l) * t_approx_P(
i,
k) + t_R(
i,
k) * t_approx_P(
i,
l));
345 t_approx_P_adjoint__dstretch(
l,
k) /= 2.;
347 t_approx_P_adjoint_log_du(
L) =
348 (t_approx_P_adjoint__dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L) +
349 t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
l,
k,
L) +
350 t_approx_P_adjoint__dstretch(
l,
k) * t_Ldiff_u(
k,
l,
L) +
351 t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
k,
l,
L)) /
355 t_levi_kirchhoff(
m) =
357 half_r * (t_diff_Rr(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
358 t_diff_Rr(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)))
362 half_l * (t_diff_Rl(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
363 t_diff_Rl(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)));
364 t_levi_kirchhoff(
m) /= 2.;
369 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
373 half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
375 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
378 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
380 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
381 t_R(
i,
l) * (t_Ldiff_u(
l,
k,
L) + t_Ldiff_u(
k,
l,
L)) / 2.;
385 t_A(
k,
l,
m) = t_diff_Rr(
i,
l,
m) * t_approx_P(
i,
k) +
386 t_diff_Rr(
i,
k,
m) * t_approx_P(
i,
l);
389 t_B(
k,
l,
m) = t_diff_Rl(
i,
l,
m) * t_approx_P(
i,
k) +
390 t_diff_Rl(
i,
k,
m) * t_approx_P(
i,
l);
393 t_approx_P_adjoint_log_du_domega(
m,
L) =
394 half_r * (t_A(
k,
l,
m) *
395 (t_Ldiff_u(
k,
l,
L) + t_Ldiff_u(
k,
l,
L)) / 2) +
396 half_l * (t_B(
k,
l,
m) *
397 (t_Ldiff_u(
k,
l,
L) + t_Ldiff_u(
k,
l,
L)) / 2);
400 t_A(
k,
l,
m) = t_diff_R(
i,
l,
m) * t_approx_P(
i,
k);
401 t_approx_P_adjoint_log_du_domega(
m,
L) =
402 t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L);
405 t_levi_kirchhoff_dstreach(
m,
L) =
407 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)))
412 (t_diff_Rl(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)));
414 t_levi_kirchhoff_dP(
m,
i,
k) =
415 half_r * t_diff_Rr(
i,
l,
m) * (t_u(
l,
k))
419 half_l * t_diff_Rl(
i,
l,
m) * (t_u(
l,
k));
421 t_levi_kirchhoff_domega(
m,
n) =
423 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
424 t_diff_diff_Rr(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)))
429 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
430 t_diff_diff_Rl(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)));
431 t_levi_kirchhoff_domega(
m,
n) /= 2.;
440 auto large_loop = [&]() {
450 "rotSelector should be large");
453 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
467 "gradApproximator not handled");
471 auto t_Ldiff_u = calculate_log_stretch();
473 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
476 t_h_u(
l,
k) = t_u(
l, o) * t_h1(o,
k);
478 t_Ldiff_h_u(
l,
k,
L) = t_Ldiff_u(
l, o,
L) * t_h1(o,
k);
494 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
495 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
497 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
498 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
501 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
502 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
504 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
505 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
507 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
508 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
513 "rotationSelector not handled");
516 constexpr double half_r = 1 / 2.;
517 constexpr double half_l = 1 / 2.;
520 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
523 t_approx_P_adjoint__dstretch(
l, o) =
524 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
525 t_approx_P_adjoint_log_du(
L) =
526 t_approx_P_adjoint__dstretch(
l, o) * t_Ldiff_u(
l, o,
L);
529 t_levi_kirchhoff(
m) =
531 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
535 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
540 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
544 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
546 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
549 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
551 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
552 t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
556 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
558 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
560 t_approx_P_adjoint_log_du_domega(
m,
L) =
561 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
565 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
568 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
569 t_approx_P_adjoint_log_du_domega(
m,
L) =
570 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
573 t_levi_kirchhoff_dstreach(
m,
L) =
575 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
579 half_l * (t_diff_Rl(
i,
l,
m) *
580 (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
582 t_levi_kirchhoff_dP(
m,
i,
k) =
584 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
588 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
590 t_levi_kirchhoff_domega(
m,
n) =
592 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
597 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
606 auto moderate_loop = [&]() {
616 "rotSelector should be large");
619 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
630 "gradApproximator not handled");
634 auto t_Ldiff_u = calculate_log_stretch();
636 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
639 t_h_u(
l,
k) = (t_symm_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
641 t_L_h(
l,
k,
L) = t_L(
l, o,
L) * t_h1(o,
k);
657 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
658 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
660 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
661 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
664 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
665 t_diff_R(
i,
j,
l) = levi_civita(
i,
j,
l);
667 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
668 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
670 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
671 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
676 "rotationSelector not handled");
679 constexpr double half_r = 1 / 2.;
680 constexpr double half_l = 1 / 2.;
683 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
686 t_approx_P_adjoint__dstretch(
l, o) =
687 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
689 t_approx_P_adjoint_log_du(
L) =
690 t_approx_P_adjoint__dstretch(
l, o) * t_L(
l, o,
L);
693 t_levi_kirchhoff(
m) =
695 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
699 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
704 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
708 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
710 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
713 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
715 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
719 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_L_h(
l,
k,
L);
721 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_L_h(
l,
k,
L);
722 t_approx_P_adjoint_log_du_domega(
m,
L) =
723 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
727 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
730 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_L_h(
l,
k,
L);
731 t_approx_P_adjoint_log_du_domega(
m,
L) =
732 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
735 t_levi_kirchhoff_dstreach(
m,
L) =
736 half_r * (t_diff_Rr(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)))
740 half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
742 t_levi_kirchhoff_dP(
m,
i,
k) =
744 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
748 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
750 t_levi_kirchhoff_domega(
m,
n) =
752 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
757 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
766 auto small_loop = [&]() {
773 "rotSelector should be small");
776 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
785 "gradApproximator not handled");
792 t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
794 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
795 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
797 t_log_u2_h1(
i,
j) = 0;
798 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
801 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
803 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
804 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
807 t_approx_P_adjoint__dstretch(
i,
j) = t_approx_P(
i,
j);
808 t_approx_P_adjoint_log_du(
L) =
809 t_approx_P_adjoint__dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
810 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
811 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
814 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
815 t_levi_kirchhoff_dstreach(
m,
L) = 0;
816 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
817 t_levi_kirchhoff_domega(
m,
n) = 0;
843 "gradApproximator not handled");
856 auto N_InLoop = getNinTheLoop();
857 auto sensee = getSkeletonSense();
858 auto nb_gauss_pts = getGaussPts().size2();
859 auto t_normal = getFTensor1NormalsAtGaussPts();
862 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
864 dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
865 dataAtPts->tractionAtPts.clear();
867 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
869 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
870 t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
882 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
888 int nb_integration_pts = getGaussPts().size2();
889 auto t_w = getFTensor0IntegrationWeight();
890 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
891 auto t_coords = getFTensor1CoordsAtGaussPts();
892 auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
900 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
901 loc_reaction_forces(
i) += (t_traction(
i)) * t_w * getMeasure();
902 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
904 loc_moment_forces(
i) +=
905 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
906 t_traction(
k) * t_w * getMeasure();
913 reactionVec[0] += loc_reaction_forces(0);
914 reactionVec[1] += loc_reaction_forces(1);
915 reactionVec[2] += loc_reaction_forces(2);
916 reactionVec[3] += loc_moment_forces(0);
917 reactionVec[4] += loc_moment_forces(1);
918 reactionVec[5] += loc_moment_forces(2);
926 int nb_integration_pts = data.
getN().size1();
927 auto v = getVolume();
928 auto t_w = getFTensor0IntegrationWeight();
929 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
930 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
931 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
932 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
933 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
934 dataAtPts->wL2DotDotAtPts.clear();
936 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
938 auto piola_scale = dataAtPts->piolaScale;
939 auto alpha_w = alphaW / piola_scale;
940 auto alpha_rho = alphaRho / piola_scale;
942 int nb_base_functions = data.
getN().size2();
946 auto get_ftensor1 = [](
auto &
v) {
951 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
953 auto t_nf = get_ftensor1(nF);
955 for (; bb != nb_dofs / 3; ++bb) {
956 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
957 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
958 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
962 for (; bb != nb_base_functions; ++bb)
976 int nb_integration_pts = getGaussPts().size2();
977 auto v = getVolume();
978 auto t_w = getFTensor0IntegrationWeight();
979 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
980 auto t_omega_grad_dot =
981 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
982 int nb_base_functions = data.
getN().size2();
988 auto get_ftensor1 = [](
auto &
v) {
994 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
997 auto t_nf = get_ftensor1(nF);
999 for (; bb != nb_dofs / 3; ++bb) {
1000 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
1001 t_nf(
k) += (
a * alphaOmega ) *
1002 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
1007 for (; bb != nb_base_functions; ++bb) {
1021 int nb_integration_pts = data.
getN().size1();
1022 auto v = getVolume();
1023 auto t_w = getFTensor0IntegrationWeight();
1025 int nb_base_functions = data.
getN().size2() / 3;
1033 auto get_ftensor1 = [](
auto &
v) {
1040 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1041 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1043 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1045 auto t_nf = get_ftensor1(nF);
1050 for (; bb != nb_dofs / 3; ++bb) {
1051 t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
1052 t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
1053 t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
1058 for (; bb != nb_base_functions; ++bb)
1068 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1070 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1072 auto t_nf = get_ftensor1(nF);
1080 for (; bb != nb_dofs / 3; ++bb) {
1081 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1086 for (; bb != nb_base_functions; ++bb)
1100 int nb_integration_pts = data.
getN().size1();
1101 auto v = getVolume();
1102 auto t_w = getFTensor0IntegrationWeight();
1104 int nb_base_functions = data.
getN().size2() / 9;
1112 auto get_ftensor0 = [](
auto &
v) {
1118 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1119 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1121 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1123 auto t_nf = get_ftensor0(nF);
1128 for (; bb != nb_dofs; ++bb) {
1129 t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
1130 t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
1134 for (; bb != nb_base_functions; ++bb) {
1144 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1146 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1148 auto t_nf = get_ftensor0(nF);
1152 t_residuum(
i,
j) = t_h(
i,
j);
1155 for (; bb != nb_dofs; ++bb) {
1156 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1160 for (; bb != nb_base_functions; ++bb) {
1174 int nb_integration_pts = data.
getN().size1();
1175 auto v = getVolume();
1176 auto t_w = getFTensor0IntegrationWeight();
1177 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1178 int nb_base_functions = data.
getN().size2() / 3;
1181 auto get_ftensor1 = [](
auto &
v) {
1186 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1188 auto t_nf = get_ftensor1(nF);
1190 for (; bb != nb_dofs / 3; ++bb) {
1191 double div_row_base = t_row_diff_base_fun(
i,
i);
1192 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1194 ++t_row_diff_base_fun;
1196 for (; bb != nb_base_functions; ++bb) {
1197 ++t_row_diff_base_fun;
1211 int nb_integration_pts = getGaussPts().size2();
1214 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1216 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1217 if (tag_length != 9) {
1219 "Number of internal stress components should be 9 but is %d",
1224 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1225 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1226 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1227 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1229 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1230 dataAtPts->internalStressAtPts.clear();
1231 auto t_internal_stress =
1232 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1235 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1236 t_internal_stress(
L) = t_const_stress(
L);
1237 ++t_internal_stress;
1249 int nb_integration_pts = getGaussPts().size2();
1252 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1254 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1255 if (tag_length != 9) {
1257 "Number of internal stress components should be 9 but is %d",
1261 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1264 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1267 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1270 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1271 dataAtPts->internalStressAtPts.clear();
1272 auto t_internal_stress =
1273 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1278 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1279 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1280 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1281 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1285 ++t_internal_stress;
1296 int nb_integration_pts = data.
getN().size1();
1297 auto v = getVolume();
1298 auto t_w = getFTensor0IntegrationWeight();
1303 auto get_ftensor2 = [](
auto &
v) {
1305 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1308 auto t_internal_stress =
1309 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1314 double scale = scalingMethodPtr->getScale(time);
1319 int nb_base_functions = data.
getN().size2();
1321 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1323 auto t_nf = get_ftensor2(nF);
1326 t_symm_stress(
i,
j) =
1327 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1330 t_residual(
L) = t_L(
i,
j,
L) * (
scale * t_symm_stress(
i,
j));
1333 for (; bb != nb_dofs / 6; ++bb) {
1334 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1338 for (; bb != nb_base_functions; ++bb)
1342 ++t_internal_stress;
1352 int nb_integration_pts = data.
getN().size1();
1353 auto v = getVolume();
1354 auto t_w = getFTensor0IntegrationWeight();
1356 auto get_ftensor2 = [](
auto &
v) {
1358 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1361 auto t_internal_stress =
1362 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1372 double scale = scalingMethodPtr->getScale(time);
1374 int nb_base_functions = data.
getN().size2();
1376 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1378 auto t_nf = get_ftensor2(nF);
1381 t_residual(
L) = t_L(M,
L) * (
scale * t_internal_stress(M));
1384 for (; bb != nb_dofs / 6; ++bb) {
1385 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1389 for (; bb != nb_base_functions; ++bb)
1393 ++t_internal_stress;
1398template <AssemblyType A>
1404 for (
auto &bc : (*bcDispPtr)) {
1406 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1409 int nb_integration_pts = OP::getGaussPts().size2();
1410 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1411 auto t_w = OP::getFTensor0IntegrationWeight();
1412 int nb_base_functions = data.
getN().size2() / 3;
1419 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1421 scale *= scalingMethodsMap.at(bc.blockName)
1424 scale *= scalingMethodsMap.at(bc.blockName)
1425 ->getScale(OP::getFEMethod()->ts_t);
1429 <<
"No scaling method found for " << bc.blockName;
1436 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1437 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1439 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1441 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1445 for (; bb != nb_base_functions; ++bb)
1457 return OP::iNtegrate(data);
1460template <AssemblyType A>
1468 double time = OP::getFEMethod()->ts_t;
1476 for (
auto &bc : (*bcRotPtr)) {
1478 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1480 int nb_integration_pts = OP::getGaussPts().size2();
1481 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1482 auto t_w = OP::getFTensor0IntegrationWeight();
1484 int nb_base_functions = data.
getN().size2() / 3;
1487 auto get_ftensor1 = [](
auto &
v) {
1500 auto get_rotation_angle = [&]() {
1501 double theta = bc.theta;
1502 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1503 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1508 auto get_rotation = [&](
auto theta) {
1510 if (bc.vals.size() == 7) {
1511 t_omega(0) = bc.vals[4];
1512 t_omega(1) = bc.vals[5];
1513 t_omega(2) = bc.vals[6];
1516 t_omega(
i) = OP::getFTensor1Normal()(
i);
1519 t_omega(
i) *= theta;
1526 auto t_R = get_rotation(get_rotation_angle());
1527 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1529 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1531 t_delta(
i) = t_center(
i) - t_coords(
i);
1533 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1535 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1537 for (; bb != nb_dofs / 3; ++bb) {
1538 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1542 for (; bb != nb_base_functions; ++bb)
1555 return OP::iNtegrate(data);
1558template <AssemblyType A>
1562 double time = OP::getFEMethod()->ts_t;
1570 for (
auto &bc : (*bcDispPtr)) {
1572 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1574 for (
auto &bd : (*brokenBaseSideDataPtr)) {
1576 auto t_approx_P = getFTensor2FromMat<3, 3>(bd.getFlux());
1577 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1578 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1579 auto t_w = OP::getFTensor0IntegrationWeight();
1587 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1588 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1591 <<
"No scaling method found for " << bc.blockName;
1595 double val =
scale * bc.val;
1598 int nb_integration_pts = OP::getGaussPts().size2();
1599 int nb_base_functions = data.
getN().size2();
1601 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1604 t_N(
i) = t_normal(
i);
1608 t_P(
i,
j) = t_N(
i) * t_N(
j);
1613 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
1617 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
1619 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1621 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1622 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
1626 for (; bb != nb_base_functions; ++bb)
1640template <AssemblyType A>
1646 double time = OP::getFEMethod()->ts_t;
1651 int row_nb_dofs = row_data.
getIndices().size();
1652 int col_nb_dofs = col_data.
getIndices().size();
1653 auto &locMat = OP::locMat;
1654 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1660 for (
auto &bc : (*bcDispPtr)) {
1662 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1664 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1665 auto t_w = OP::getFTensor0IntegrationWeight();
1671 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1672 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1675 <<
"No scaling method found for " << bc.blockName;
1678 int nb_integration_pts = OP::getGaussPts().size2();
1679 int row_nb_dofs = row_data.
getIndices().size();
1680 int col_nb_dofs = col_data.
getIndices().size();
1681 int nb_base_functions = row_data.
getN().size2();
1686 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1689 t_N(
i) = t_normal(
i);
1693 t_P(
i,
j) = t_N(
i) * t_N(
j);
1696 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1699 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1700 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1703 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1704 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1711 for (; rr != nb_base_functions; ++rr)
1718 locMat *= OP::getMeasure();
1725template <AssemblyType A>
1731 double time = OP::getFEMethod()->ts_t;
1736 int row_nb_dofs = row_data.
getIndices().size();
1737 int col_nb_dofs = col_data.
getIndices().size();
1738 auto &locMat = OP::locMat;
1739 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1745 for (
auto &bc : (*bcDispPtr)) {
1747 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1749 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1750 auto t_w = OP::getFTensor0IntegrationWeight();
1759 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1760 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1763 <<
"No scaling method found for " << bc.blockName;
1766 int nb_integration_pts = OP::getGaussPts().size2();
1767 int nb_base_functions = row_data.
getN().size2();
1770 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1773 t_N(
i) = t_normal(
i);
1777 t_P(
i,
j) = t_N(
i) * t_N(
j);
1782 t_d_res(
i,
j) = t_Q(
i,
j);
1785 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1786 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1789 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1791 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1798 for (; rr != nb_base_functions; ++rr)
1805 locMat *= OP::getMeasure();
1812 return OP::iNtegrate(data);
1817 return OP::iNtegrate(row_data, col_data);
1822 return OP::iNtegrate(row_data, col_data);
1825template <AssemblyType A>
1829 double time = OP::getFEMethod()->ts_t;
1837 for (
auto &bc : (*bcDispPtr)) {
1839 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1844 auto [block_name, v_analytical_expr] =
1849 int nb_integration_pts = OP::getGaussPts().size2();
1850 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1851 auto t_w = OP::getFTensor0IntegrationWeight();
1852 int nb_base_functions = data.
getN().size2() / 3;
1854 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1860 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1862 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1863 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1866 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1868 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1872 for (; bb != nb_base_functions; ++bb)
1886 return OP::iNtegrate(data);
1895 int nb_integration_pts = getGaussPts().size2();
1896 int nb_base_functions = data.
getN().size2();
1898 double time = getFEMethod()->ts_t;
1904 if (this->locF.size() != nb_dofs)
1906 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1909 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1912 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1914 auto t_w = getFTensor0IntegrationWeight();
1915 auto t_coords = getFTensor1CoordsAtGaussPts();
1917 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1919 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1921 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1922 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1924 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1925 t_f(
i) -= (time_scale * t_w * t_row_base * tau) * (t_val(
i) *
scale);
1930 for (; rr != nb_base_functions; ++rr)
1935 this->locF *= getMeasure();
1941 for (
auto &bc : *(bcData)) {
1942 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1944 double time_scale = 1;
1945 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1946 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1952 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1956 return -y * (y - 1) / 0.25;
1958 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1961 bc, [](
double,
double,
double) {
return 1; }, time_scale);
1975 int nb_integration_pts = getGaussPts().size2();
1976 int nb_base_functions = data.
getN().size2();
1978 double time = getFEMethod()->ts_t;
1984 if (this->locF.size() != nb_dofs)
1986 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1989 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1994 auto t_w = getFTensor0IntegrationWeight();
1995 auto t_coords = getFTensor1CoordsAtGaussPts();
1996 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
1997 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
1999 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2001 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2003 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2015 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2018 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2019 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2020 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2022 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2024 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2026 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2028 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2029 t_f(
i) += t_row_base * t_val(
i);
2034 for (; rr != nb_base_functions; ++rr)
2049 for (
auto &bc : *(bcData)) {
2050 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2052 double time_scale = 1;
2053 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2054 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2060 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2067template <AssemblyType A>
2078 double time = OP::getFEMethod()->ts_t;
2083 int nb_base_functions = row_data.
getN().size2();
2084 int row_nb_dofs = row_data.
getIndices().size();
2085 int col_nb_dofs = col_data.
getIndices().size();
2086 int nb_integration_pts = OP::getGaussPts().size2();
2087 auto &locMat = OP::locMat;
2088 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2091auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2096 auto t_w = OP::getFTensor0IntegrationWeight();
2097 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2098 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2099 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2101 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2104 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2114 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2115 auto t_val = time_scale * t_w * tau * val;
2118 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2119 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2122 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2124 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2125 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2126 t_kd(
j,
l) * t_diff_col_base(N0)
2130 (FTensor::levi_civita<double>(
i,
j,
k) *
2131 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2132 t_kd(
k,
l) * t_diff_col_base(N1);
2134 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2141 for (; rr != nb_base_functions; ++rr)
2157 for (
auto &bc : *(bcData)) {
2158 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2160 double time_scale = 1;
2161 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2162 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2168 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2179 return OP::iNtegrate(row_data, col_data);
2189 int nb_integration_pts = getGaussPts().size2();
2190 int nb_base_functions = data.
getN().size2();
2192 double time = getFEMethod()->ts_t;
2198 if (this->locF.size() != nb_dofs)
2200 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2205 for (
auto &bc : *(bcData)) {
2206 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2210 auto [block_name, v_analytical_expr] =
2213 getFTensor1FromMat<3>(v_analytical_expr);
2215 auto t_w = getFTensor0IntegrationWeight();
2216 auto t_coords = getFTensor1CoordsAtGaussPts();
2218 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2220 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2222 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2224 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2225 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2230 for (; rr != nb_base_functions; ++rr)
2236 this->locF *= getMeasure();
2245 int nb_integration_pts = row_data.
getN().size1();
2246 int row_nb_dofs = row_data.
getIndices().size();
2247 int col_nb_dofs = col_data.
getIndices().size();
2248 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2250 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2253 auto v = getVolume();
2254 auto t_w = getFTensor0IntegrationWeight();
2255 int row_nb_base_functions = row_data.
getN().size2();
2257 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2260 for (; rr != row_nb_dofs / 3; ++rr) {
2262 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2263 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2264 double div_col_base = t_col_diff_base_fun(
i,
i);
2265 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2267 ++t_col_diff_base_fun;
2271 for (; rr != row_nb_base_functions; ++rr)
2282 if (alphaW < std::numeric_limits<double>::epsilon() &&
2283 alphaRho < std::numeric_limits<double>::epsilon())
2286 const int nb_integration_pts = row_data.
getN().size1();
2287 const int row_nb_dofs = row_data.
getIndices().size();
2288 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2290 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2296 auto v = getVolume();
2297 auto t_w = getFTensor0IntegrationWeight();
2299 auto piola_scale = dataAtPts->piolaScale;
2300 auto alpha_w = alphaW / piola_scale;
2301 auto alpha_rho = alphaRho / piola_scale;
2303 int row_nb_base_functions = row_data.
getN().size2();
2306 double ts_scale = alpha_w * getTSa();
2307 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2308 ts_scale += alpha_rho * getTSaa();
2310 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2311 double a =
v * t_w * ts_scale;
2314 for (; rr != row_nb_dofs / 3; ++rr) {
2317 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2318 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2319 const double b =
a * t_row_base_fun * t_col_base_fun;
2328 for (; rr != row_nb_base_functions; ++rr)
2347 int nb_integration_pts = row_data.
getN().size1();
2348 int row_nb_dofs = row_data.
getIndices().size();
2349 int col_nb_dofs = col_data.
getIndices().size();
2350 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2353 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2355 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2357 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2359 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2361 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2363 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2366 auto v = getVolume();
2367 auto t_w = getFTensor0IntegrationWeight();
2369 int row_nb_base_functions = row_data.
getN().size2();
2372 auto t_approx_P_adjoint_log_du_dP =
2373 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2375 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2378 for (; rr != row_nb_dofs / 6; ++rr) {
2381 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2383 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2385 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2393 for (; rr != row_nb_base_functions; ++rr)
2396 ++t_approx_P_adjoint_log_du_dP;
2412 int nb_integration_pts = row_data.
getN().size1();
2413 int row_nb_dofs = row_data.
getIndices().size();
2414 int col_nb_dofs = col_data.
getIndices().size();
2415 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2417 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2421 auto v = getVolume();
2422 auto t_w = getFTensor0IntegrationWeight();
2425 int row_nb_base_functions = row_data.
getN().size2();
2427 auto t_approx_P_adjoint_log_du_dP =
2428 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2430 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2433 for (; rr != row_nb_dofs / 6; ++rr) {
2434 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2435 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2436 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2438 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2445 for (; rr != row_nb_base_functions; ++rr)
2448 ++t_approx_P_adjoint_log_du_dP;
2460 int row_nb_dofs = row_data.
getIndices().size();
2461 int col_nb_dofs = col_data.
getIndices().size();
2462 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2465 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2467 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2469 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2471 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2473 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2475 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2485 auto v = getVolume();
2486 auto t_w = getFTensor0IntegrationWeight();
2487 auto t_approx_P_adjoint_log_du_domega =
2488 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2490 int row_nb_base_functions = row_data.
getN().size2();
2493 int nb_integration_pts = row_data.
getN().size1();
2495 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2499 for (; rr != row_nb_dofs / 6; ++rr) {
2501 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2502 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2503 double v =
a * t_row_base_fun * t_col_base_fun;
2504 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2511 for (; rr != row_nb_base_functions; ++rr)
2515 ++t_approx_P_adjoint_log_du_domega;
2524 int nb_integration_pts = getGaussPts().size2();
2525 int row_nb_dofs = row_data.
getIndices().size();
2526 int col_nb_dofs = col_data.
getIndices().size();
2527 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2531 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2532 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2534 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2535 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2537 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2538 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2546 auto v = getVolume();
2547 auto t_w = getFTensor0IntegrationWeight();
2548 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2549 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2550 int row_nb_base_functions = row_data.
getN().size2();
2552 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2555 for (; rr != row_nb_dofs / 3; ++rr) {
2556 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2557 const double b =
a * t_row_base_fun;
2559 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2560 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2566 for (; rr != row_nb_base_functions; ++rr) {
2570 ++t_levi_kirchhoff_du;
2586 int nb_integration_pts = getGaussPts().size2();
2587 int row_nb_dofs = row_data.
getIndices().size();
2588 int col_nb_dofs = col_data.
getIndices().size();
2589 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2592 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2594 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2596 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2601 auto v = getVolume();
2602 auto t_w = getFTensor0IntegrationWeight();
2604 int row_nb_base_functions = row_data.
getN().size2();
2606 auto t_levi_kirchhoff_dP =
2607 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2609 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2612 for (; rr != row_nb_dofs / 3; ++rr) {
2613 double b =
a * t_row_base_fun;
2615 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2616 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2617 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2623 for (; rr != row_nb_base_functions; ++rr) {
2628 ++t_levi_kirchhoff_dP;
2636 int nb_integration_pts = getGaussPts().size2();
2637 int row_nb_dofs = row_data.
getIndices().size();
2638 int col_nb_dofs = col_data.
getIndices().size();
2640 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2642 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2649 auto v = getVolume();
2650 auto t_w = getFTensor0IntegrationWeight();
2651 auto t_levi_kirchoff_dP =
2652 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2654 int row_nb_base_functions = row_data.
getN().size2();
2657 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2660 for (; rr != row_nb_dofs / 3; ++rr) {
2661 double b =
a * t_row_base_fun;
2662 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2663 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2664 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2665 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2672 for (; rr != row_nb_base_functions; ++rr) {
2676 ++t_levi_kirchoff_dP;
2684 int nb_integration_pts = getGaussPts().size2();
2685 int row_nb_dofs = row_data.
getIndices().size();
2686 int col_nb_dofs = col_data.
getIndices().size();
2687 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2690 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2692 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2694 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2707 auto v = getVolume();
2708 auto ts_a = getTSa();
2709 auto t_w = getFTensor0IntegrationWeight();
2710 auto t_levi_kirchhoff_domega =
2711 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2712 int row_nb_base_functions = row_data.
getN().size2();
2717 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2719 double c = (
a * alphaOmega) * (ts_a );
2722 for (; rr != row_nb_dofs / 3; ++rr) {
2723 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2724 const double b =
a * t_row_base_fun;
2727 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2728 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2729 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2737 for (; rr != row_nb_base_functions; ++rr) {
2742 ++t_levi_kirchhoff_domega;
2750 if (dataAtPts->matInvD.size1() ==
size_symm &&
2751 dataAtPts->matInvD.size2() ==
size_symm) {
2764 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2767 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2769 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2771 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2776 int nb_integration_pts = getGaussPts().size2();
2777 int row_nb_dofs = row_data.
getIndices().size();
2778 int col_nb_dofs = col_data.
getIndices().size();
2780 auto v = getVolume();
2781 auto t_w = getFTensor0IntegrationWeight();
2782 int row_nb_base_functions = row_data.
getN().size2() / 3;
2789 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2790 &*dataAtPts->matInvD.data().begin());
2793 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2797 for (; rr != row_nb_dofs / 3; ++rr) {
2799 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2800 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2801 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2809 for (; rr != row_nb_base_functions; ++rr)
2822 if (dataAtPts->matInvD.size1() ==
size_symm &&
2823 dataAtPts->matInvD.size2() ==
size_symm) {
2837 int nb_integration_pts = getGaussPts().size2();
2838 int row_nb_dofs = row_data.
getIndices().size();
2839 int col_nb_dofs = col_data.
getIndices().size();
2841 auto v = getVolume();
2842 auto t_w = getFTensor0IntegrationWeight();
2843 int row_nb_base_functions = row_data.
getN().size2() / 9;
2850 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2851 &*dataAtPts->matInvD.data().begin());
2854 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2858 for (; rr != row_nb_dofs; ++rr) {
2860 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2862 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2869 for (; rr != row_nb_base_functions; ++rr)
2880 if (dataAtPts->matInvD.size1() ==
size_symm &&
2881 dataAtPts->matInvD.size2() ==
size_symm) {
2895 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2898 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
2903 int nb_integration_pts = getGaussPts().size2();
2904 int row_nb_dofs = row_data.
getIndices().size();
2905 int col_nb_dofs = col_data.
getIndices().size();
2907 auto v = getVolume();
2908 auto t_w = getFTensor0IntegrationWeight();
2909 int row_nb_base_functions = row_data.
getN().size2() / 9;
2918 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2919 &*dataAtPts->matInvD.data().begin());
2922 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2925 auto t_m = get_ftensor1(
K, 0, 0);
2928 for (; rr != row_nb_dofs; ++rr) {
2930 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2931 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2939 for (; rr != row_nb_base_functions; ++rr)
2958 int nb_integration_pts = row_data.
getN().size1();
2959 int row_nb_dofs = row_data.
getIndices().size();
2960 int col_nb_dofs = col_data.
getIndices().size();
2962 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2965 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2967 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2969 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2974 auto v = getVolume();
2975 auto t_w = getFTensor0IntegrationWeight();
2976 int row_nb_base_functions = row_data.
getN().size2() / 3;
2979 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2981 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2985 for (; rr != row_nb_dofs / 3; ++rr) {
2988 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2991 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2992 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2993 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
3001 for (; rr != row_nb_base_functions; ++rr)
3021 int nb_integration_pts = row_data.
getN().size1();
3022 int row_nb_dofs = row_data.
getIndices().size();
3023 int col_nb_dofs = col_data.
getIndices().size();
3025 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3027 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3030 auto v = getVolume();
3031 auto t_w = getFTensor0IntegrationWeight();
3032 int row_nb_base_functions = row_data.
getN().size2() / 9;
3035 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3036 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3040 for (; rr != row_nb_dofs; ++rr) {
3043 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3046 auto t_m = get_ftensor2(
K, rr, 0);
3047 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3048 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3056 for (; rr != row_nb_base_functions; ++rr)
3069 if (tagSense != getSkeletonSense())
3072 auto get_tag = [&](
auto name) {
3073 auto &mob = getPtrFE()->mField.get_moab();
3079 auto get_tag_value = [&](
auto &&tag,
int dim) {
3080 auto &mob = getPtrFE()->mField.get_moab();
3081 auto face = getSidePtrFE()->getFEEntityHandle();
3082 std::vector<double> value(dim);
3083 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3087 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3088 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3090 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3091 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3096 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3097 Tag th_detF = create_tag(
"detF", 1);
3098 Tag th_traction = create_tag(
"traction", 3);
3099 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3101 Tag th_energy = create_tag(
"Energy", 1);
3103 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3104 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3105 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3107 auto t_normal = getFTensor1NormalsAtGaussPts();
3108 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3120 if (dataAtPts->energyAtPts.size() == 0) {
3122 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3123 dataAtPts->energyAtPts.clear();
3132 auto set_float_precision = [](
const double x) {
3133 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3140 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3142 v = set_float_precision(
v);
3143 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3150 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3153 for (
auto &
a :
v.data())
3154 a = set_float_precision(
a);
3155 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3156 &*
v.data().begin());
3164 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3166 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3168 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3170 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3172 t_m(
i,
j) = t_d(
i,
j);
3173 for (
auto &
v :
m.data())
3174 v = set_float_precision(
v);
3175 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3176 &*
m.data().begin());
3180 const auto nb_gauss_pts = getGaussPts().size2();
3181 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3184 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3186 t_traction(
i) *= tagSense;
3187 CHKERR save_vec_tag(th_traction, t_traction, gg);
3189 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3190 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3191 CHKERR save_scal_tag(th_energy, t_energy, gg);
3195 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3196 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3197 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3206 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3207 std::vector<FieldSpace> spaces, std::string geom_field_name,
3208 boost::shared_ptr<Range> crack_front_edges_ptr) {
3211 constexpr bool scale_l2 =
false;
3215 "Scale L2 Ainsworth Legendre base is not implemented");
3224 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3225 std::vector<FieldSpace> spaces, std::string geom_field_name,
3226 boost::shared_ptr<Range> crack_front_edges_ptr) {
3229 constexpr bool scale_l2 =
false;
3233 "Scale L2 Ainsworth Legendre base is not implemented");
3242 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3243 std::vector<FieldSpace> spaces, std::string geom_field_name,
3244 boost::shared_ptr<Range> crack_front_edges_ptr,
3245 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3246 boost::shared_ptr<MatrixDouble> inv_jac) {
3250 auto jac = boost::make_shared<MatrixDouble>();
3251 auto det = boost::make_shared<VectorDouble>();
3253 geom_field_name, jac));
3259 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3261 if (scale_l2_ainsworth_legendre_base) {
3269 boost::shared_ptr<MatrixDouble> jac,
3270 boost::shared_ptr<Range> edges_ptr)
3279 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3282 return OP::doWork(side, type, data);
3287 boost::shared_ptr<Range> edgesPtr;
3290 auto jac = boost::make_shared<MatrixDouble>();
3291 auto det = boost::make_shared<VectorDouble>();
3293 geom_field_name, jac,
3295 : boost::make_shared<
Range>()));
3365 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3366 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3367 if (getNinTheLoop() == 0) {
3368 dataAtPts->faceMaterialForceAtPts.clear();
3369 dataAtPts->normalPressureAtPts.clear();
3371 auto loop_size = getLoopSize();
3372 if (loop_size == 1) {
3373 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3374 auto pstatus = numebered_fe_ptr->getPStatus();
3375 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3382 auto t_normal = getFTensor1NormalsAtGaussPts();
3383 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3384 dataAtPts->faceMaterialForceAtPts);
3387 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3391 getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
3392 auto t_grad_u_gamma =
3393 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3395 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3396 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
3418 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3419 t_N(
I) = t_normal(
I);
3422 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3424 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3427 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3430 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3432 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3438 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3441 t_N(
I) = t_normal(
I);
3446 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3449 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3450 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3453 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3456 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3459 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3467 "Grffith energy release "
3468 "selector not implemented");
3472 auto side_fe_ptr = getSidePtrFE();
3473 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3474 auto pstatus = side_fe_mi_ptr->getPStatus();
3476 auto owner = side_fe_mi_ptr->getOwnerProc();
3478 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3479 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3480 << getNinTheLoop() <<
" loop size " << getLoopSize();
3492 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3493 auto pstatus = fe_mi_ptr->getPStatus();
3495 auto owner = fe_mi_ptr->getOwnerProc();
3497 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3498 <<
" " << getPtrFE()->mField.get_comm_rank();
3506 double face_pressure = 0.;
3507 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3508 dataAtPts->faceMaterialForceAtPts);
3511 auto t_w = getFTensor0IntegrationWeight();
3512 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3513 t_face_T(
I) += t_w * t_T(
I);
3514 face_pressure += t_w * t_p;
3519 t_face_T(
I) *= getMeasure();
3520 face_pressure *= getMeasure();
3522 auto get_tag = [&](
auto name,
auto dim) {
3523 auto &moab = getPtrFE()->mField.get_moab();
3525 double def_val[] = {0., 0., 0.};
3526 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3527 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3532 auto set_tag = [&](
auto &&tag,
auto ptr) {
3533 auto &moab = getPtrFE()->mField.get_moab();
3534 auto face = getPtrFE()->getFEEntityHandle();
3535 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3538 set_tag(
get_tag(
"MaterialForce", 3), &t_face_T(0));
3539 set_tag(
get_tag(
"FacePressure", 1), &face_pressure);
3545template <
typename OP_PTR>
3547 const std::string block_name) {
3549 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
3551 auto ts_time = op_ptr->getTStime();
3552 auto ts_time_step = op_ptr->getTStimeStep();
3558 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
3559 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
3561 auto v_analytical_expr =
3563 m_ref_coords, m_ref_normals, block_name);
3566 if (v_analytical_expr.size2() != nb_gauss_pts)
3568 "Wrong number of integration pts");
3571 return std::make_tuple(block_name, v_analytical_expr);
3575#ifdef ENABLE_PYTHON_BINDING
3577 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3582 auto main_module = bp::import(
"__main__");
3583 mainNamespace = main_module.attr(
"__dict__");
3584 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3586 analyticalDispFun = mainNamespace[
"analytical_disp"];
3587 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3588 analyticalExternalStrainFun = mainNamespace[
"analytical_external_strain"];
3590 }
catch (bp::error_already_set
const &) {
3600 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3601 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3602 const std::string &block_name, np::ndarray &analytical_expr
3609 analytical_expr = bp::extract<np::ndarray>(
3610 analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3612 }
catch (bp::error_already_set
const &) {
3622 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3623 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3624 const std::string& block_name, np::ndarray &analytical_expr
3631 analytical_expr = bp::extract<np::ndarray>(
3632 analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3634 }
catch (bp::error_already_set
const &) {
3642 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalExternalStrain(
3644 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3645 const std::string& block_name, np::ndarray &analytical_expr
3652 analytical_expr = bp::extract<np::ndarray>(
3653 analyticalExternalStrainFun(delta_t,
t, x, y, z, block_name));
3655 }
catch (bp::error_already_set
const &) {
3663boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3665inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3667 auto dtype = np::dtype::get_builtin<double>();
3668 auto size = bp::make_tuple(nb_gauss_pts);
3669 auto stride = bp::make_tuple(3 *
sizeof(
double));
3670 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3679 const std::string block_name) {
3681#ifdef ENABLE_PYTHON_BINDING
3682 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3687 bp::list python_coords;
3688 bp::list python_normals;
3690 for (
int idx = 0; idx < 3; ++idx) {
3691 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3692 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3695 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3696 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3698 std::regex reg_disp_name(disp_block_name);
3699 std::regex reg_traction_name(traction_block_name);
3701 np::ndarray np_analytical_expr = np::empty(
3702 bp::make_tuple(nb_gauss_pts, 3), np::dtype::get_builtin<double>());
3704 if (std::regex_match(block_name, reg_disp_name)) {
3706 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3707 bp::extract<np::ndarray>(python_coords[1]),
3708 bp::extract<np::ndarray>(python_coords[2]),
3709 bp::extract<np::ndarray>(python_normals[0]),
3710 bp::extract<np::ndarray>(python_normals[1]),
3711 bp::extract<np::ndarray>(python_normals[2]),
3712 block_name, np_analytical_expr),
3713 "Failed analytical_disp() python call");
3714 }
else if (std::regex_match(block_name, reg_traction_name)) {
3716 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3717 bp::extract<np::ndarray>(python_coords[1]),
3718 bp::extract<np::ndarray>(python_coords[2]),
3719 bp::extract<np::ndarray>(python_normals[0]),
3720 bp::extract<np::ndarray>(python_normals[1]),
3721 bp::extract<np::ndarray>(python_normals[2]),
3722 block_name, np_analytical_expr),
3723 "Failed analytical_traction() python call");
3726 "Unknown analytical block");
3729 double *analytical_expr_val_ptr =
3730 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3733 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3734 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3735 for (
int idx = 0; idx < 3; ++idx)
3736 v_analytical_expr(idx, gg) =
3737 *(analytical_expr_val_ptr + (3 * gg + idx));
3739 return v_analytical_expr;
3750 const std::string block_name) {
3752#ifdef ENABLE_PYTHON_BINDING
3753 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3757 bp::list python_coords;
3759 for (
int idx = 0; idx < 3; ++idx) {
3760 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3763 auto externalstrain_block_name =
"(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
3765 std::regex reg_externalstrain_name(externalstrain_block_name);
3767 np::ndarray np_analytical_expr = np::empty(
3768 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
3770 if (std::regex_match(block_name, reg_externalstrain_name)) {
3771 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalExternalStrain(
3772 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3773 bp::extract<np::ndarray>(python_coords[1]),
3774 bp::extract<np::ndarray>(python_coords[2]), block_name,
3775 np_analytical_expr),
3776 "Failed analytical_external_strain() python call");
3779 "Unknown analytical block");
3782 double *analytical_expr_val_ptr =
3783 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3786 v_analytical_expr.resize(nb_gauss_pts,
false);
3787 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg)
3788 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
3790 return v_analytical_expr;
3799 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3800 boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
3801 boost::shared_ptr<Range> ents_ptr)
3802 :
OP(broken_base_side_data, ents_ptr) {
3803 this->sourceVec = vec;
3804 this->betaCoeff = beta_coeff;
3808 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3809 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3810 :
OP(broken_base_side_data, ents_ptr) {
3811 this->sourceVec = boost::shared_ptr<MatrixDouble>();
3812 this->betaCoeff = beta_coeff;
3821 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3826 if (!brokenBaseSideData) {
3831 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
3839 OP::nbIntegrationPts = OP::getGaussPts().size2();
3841 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
3843 OP::locF.resize(OP::nbRows,
false);
3846 CHKERR this->iNtegrate(row_data);
3848 CHKERR this->aSsemble(row_data);
3852 switch (OP::opType) {
3854 for (
auto &bd : *brokenBaseSideData) {
3856 boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
3857 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
3858 this->sourceVec.reset();
3863 (std::string(
"wrong op type ") +
3864 OpBaseDerivativesBase::OpTypeNames[OP::opType])
3872 const std::string row_field,
3873 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3874 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3875 :
OP(row_field, boost::shared_ptr<
MatrixDouble>(), beta_coeff, ents_ptr),
3876 brokenBaseSideDataPtr(broken_base_side_data) {
3877 this->betaCoeff = beta_coeff;
3883 for (
auto &bd : (*brokenBaseSideDataPtr)) {
3888 if (this->sourceVec->size1() !=
SPACE_DIM) {
3890 "Inconsistent size of the source vector");
3892 if (this->sourceVec->size2() != OP::getGaussPts().size2()) {
3894 "Inconsistent size of the source vector");
3898 CHKERR OP::iNtegrate(data);
3900 this->sourceVec.reset();
3906 std::string row_field,
3907 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3908 ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
3909 boost::shared_ptr<Range> ents_ptr)
3910 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
3912 this->betaCoeff = beta;
3917 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3918 ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
3919 :
OP(broken_base_side_data, ents_ptr) {
3921 this->betaCoeff = beta;
3930 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3935 if (!brokenBaseSideData) {
3940 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
3941 EntityType col_type,
3946 auto check_if_assemble_transpose = [&] {
3948 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
3952 }
else if (OP::assembleTranspose) {
3958 OP::rowSide = row_side;
3959 OP::rowType = row_type;
3960 OP::colSide = col_side;
3961 OP::colType = col_type;
3963 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
3965 CHKERR this->iNtegrate(row_data, col_data);
3966 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
3970 switch (OP::opType) {
3973 for (
auto &bd : *brokenBaseSideData) {
3976 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
3978 "base functions not set");
3982 OP::nbRows = bd.getData().getIndices().size();
3985 OP::nbIntegrationPts = OP::getGaussPts().size2();
3986 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
3994 bd.getSide(), bd.getSide(),
3997 bd.getType(), bd.getType(),
4000 bd.getData(), bd.getData()
4009 (std::string(
"wrong op type ") +
4010 OpBaseDerivativesBase::OpTypeNames[OP::opType])
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
Tensor1< T, Tensor_Dim > normalize()
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
VectorDouble analytical_externalstrain_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, const std::string block_name)
auto sort_eigen_vals(A &eig, B &eigen_vec)
static Tag get_tag(moab::Interface &moab, std::string tag_name, int size)
std::tuple< std::string, MatrixDouble > getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
static constexpr auto size_symm
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, MatrixDouble &m_ref_normals, const std::string block_name)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr IntegrationType I
constexpr double t
plate stiffness
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static double dynamicTime
static PetscBool l2UserBaseScale
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool dynamicRelaxation
static enum SymmetrySelector symmetrySelector
static boost::function< double(const double)> f
static PetscBool setSingularity
static boost::function< double(const double)> d_f
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_f
static auto diffExp(A &&t_w_vee, B &&theta)
static auto exp(A &&t_w_vee, B &&theta)
Add operators pushing bases from local to physical configuration.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
Get DOF values on entity.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Operator for inverting matrices at integration points.
Scale base functions by inverses of measure of element.
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseBrokenBase(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseTimesBrokenDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBrokenBaseTimesHybridDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Caluclate face material force and normal pressure at gauss points.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpHybridBaseTimesBrokenDisp(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenBaseSideDataPtr
OpHyrbridBaseBrokenBase(std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
