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) {
958 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
960 auto t_nf = get_ftensor1(nF);
962 for (; bb != nb_dofs / 3; ++bb) {
963 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
964 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
965 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
969 for (; bb != nb_base_functions; ++bb)
980 int nb_integration_pts = getGaussPts().size2();
981 auto v = getVolume();
982 auto t_w = getFTensor0IntegrationWeight();
983 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
984 auto t_omega_grad_dot =
985 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
986 int nb_base_functions = data.
getN().size2();
992 auto get_ftensor1 = [](
auto &
v) {
998 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1001 auto t_nf = get_ftensor1(nF);
1003 for (; bb != nb_dofs / 3; ++bb) {
1004 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
1005 t_nf(
k) += (
a * alphaOmega ) *
1006 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
1011 for (; bb != nb_base_functions; ++bb) {
1025 int nb_integration_pts = data.
getN().size1();
1026 auto v = getVolume();
1027 auto t_w = getFTensor0IntegrationWeight();
1029 int nb_base_functions = data.
getN().size2() / 3;
1037 auto get_ftensor1 = [](
auto &
v) {
1044 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1045 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1047 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1049 auto t_nf = get_ftensor1(nF);
1054 for (; bb != nb_dofs / 3; ++bb) {
1055 t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
1056 t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
1057 t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
1062 for (; bb != nb_base_functions; ++bb)
1072 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1074 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1076 auto t_nf = get_ftensor1(nF);
1084 for (; bb != nb_dofs / 3; ++bb) {
1085 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1090 for (; bb != nb_base_functions; ++bb)
1104 int nb_integration_pts = data.
getN().size1();
1105 auto v = getVolume();
1106 auto t_w = getFTensor0IntegrationWeight();
1108 int nb_base_functions = data.
getN().size2() / 9;
1116 auto get_ftensor0 = [](
auto &
v) {
1122 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1123 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1125 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1127 auto t_nf = get_ftensor0(nF);
1132 for (; bb != nb_dofs; ++bb) {
1133 t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
1134 t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
1138 for (; bb != nb_base_functions; ++bb) {
1148 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1150 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1152 auto t_nf = get_ftensor0(nF);
1156 t_residuum(
i,
j) = t_h(
i,
j);
1159 for (; bb != nb_dofs; ++bb) {
1160 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1164 for (; bb != nb_base_functions; ++bb) {
1178 int nb_integration_pts = data.
getN().size1();
1179 auto v = getVolume();
1180 auto t_w = getFTensor0IntegrationWeight();
1181 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1182 int nb_base_functions = data.
getN().size2() / 3;
1185 auto get_ftensor1 = [](
auto &
v) {
1190 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1192 auto t_nf = get_ftensor1(nF);
1194 for (; bb != nb_dofs / 3; ++bb) {
1195 double div_row_base = t_row_diff_base_fun(
i,
i);
1196 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1198 ++t_row_diff_base_fun;
1200 for (; bb != nb_base_functions; ++bb) {
1201 ++t_row_diff_base_fun;
1215 int nb_integration_pts = getGaussPts().size2();
1218 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1220 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1221 if (tag_length != 9) {
1223 "Number of internal stress components should be 9 but is %d",
1228 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1229 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1230 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1231 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1233 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1234 dataAtPts->internalStressAtPts.clear();
1235 auto t_internal_stress =
1236 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1239 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1240 t_internal_stress(
L) = t_const_stress(
L);
1241 ++t_internal_stress;
1253 int nb_integration_pts = getGaussPts().size2();
1256 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1258 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1259 if (tag_length != 9) {
1261 "Number of internal stress components should be 9 but is %d",
1265 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1268 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1271 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1274 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1275 dataAtPts->internalStressAtPts.clear();
1276 auto t_internal_stress =
1277 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1282 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1283 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1284 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1285 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1289 ++t_internal_stress;
1300 int nb_integration_pts = data.
getN().size1();
1301 auto v = getVolume();
1302 auto t_w = getFTensor0IntegrationWeight();
1307 auto get_ftensor2 = [](
auto &
v) {
1309 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1312 auto t_internal_stress =
1313 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1318 double scale = scalingMethodPtr->getScale(time);
1323 int nb_base_functions = data.
getN().size2();
1325 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1327 auto t_nf = get_ftensor2(nF);
1330 t_symm_stress(
i,
j) =
1331 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1334 t_residual(
L) = t_L(
i,
j,
L) * (
scale * t_symm_stress(
i,
j));
1337 for (; bb != nb_dofs / 6; ++bb) {
1338 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1342 for (; bb != nb_base_functions; ++bb)
1346 ++t_internal_stress;
1356 int nb_integration_pts = data.
getN().size1();
1357 auto v = getVolume();
1358 auto t_w = getFTensor0IntegrationWeight();
1360 auto get_ftensor2 = [](
auto &
v) {
1362 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1365 auto t_internal_stress =
1366 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1376 double scale = scalingMethodPtr->getScale(time);
1378 int nb_base_functions = data.
getN().size2();
1380 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1382 auto t_nf = get_ftensor2(nF);
1385 t_residual(
L) = t_L(M,
L) * (
scale * t_internal_stress(M));
1388 for (; bb != nb_dofs / 6; ++bb) {
1389 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1393 for (; bb != nb_base_functions; ++bb)
1397 ++t_internal_stress;
1402template <AssemblyType A>
1408 for (
auto &bc : (*bcDispPtr)) {
1410 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1413 int nb_integration_pts = OP::getGaussPts().size2();
1414 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1415 auto t_w = OP::getFTensor0IntegrationWeight();
1416 int nb_base_functions = data.
getN().size2() / 3;
1423 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1425 scale *= scalingMethodsMap.at(bc.blockName)
1428 scale *= scalingMethodsMap.at(bc.blockName)
1429 ->getScale(OP::getFEMethod()->ts_t);
1433 <<
"No scaling method found for " << bc.blockName;
1440 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1441 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1443 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1445 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1449 for (; bb != nb_base_functions; ++bb)
1461 return OP::iNtegrate(data);
1464template <AssemblyType A>
1472 double time = OP::getFEMethod()->ts_t;
1480 for (
auto &bc : (*bcRotPtr)) {
1482 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1484 int nb_integration_pts = OP::getGaussPts().size2();
1485 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1486 auto t_w = OP::getFTensor0IntegrationWeight();
1488 int nb_base_functions = data.
getN().size2() / 3;
1491 auto get_ftensor1 = [](
auto &
v) {
1504 auto get_rotation_angle = [&]() {
1505 double theta = bc.theta;
1506 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1507 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1512 auto get_rotation = [&](
auto theta) {
1514 if (bc.vals.size() == 7) {
1515 t_omega(0) = bc.vals[4];
1516 t_omega(1) = bc.vals[5];
1517 t_omega(2) = bc.vals[6];
1520 t_omega(
i) = OP::getFTensor1Normal()(
i);
1523 t_omega(
i) *= theta;
1530 auto t_R = get_rotation(get_rotation_angle());
1531 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1533 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1535 t_delta(
i) = t_center(
i) - t_coords(
i);
1537 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1539 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1541 for (; bb != nb_dofs / 3; ++bb) {
1542 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1546 for (; bb != nb_base_functions; ++bb)
1559 return OP::iNtegrate(data);
1562template <AssemblyType A>
1566 double time = OP::getFEMethod()->ts_t;
1574 for (
auto &bc : (*bcDispPtr)) {
1576 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1578 for (
auto &bd : (*brokenBaseSideDataPtr)) {
1580 auto t_approx_P = getFTensor2FromMat<3, 3>(bd.getFlux());
1581 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1582 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1583 auto t_w = OP::getFTensor0IntegrationWeight();
1591 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1592 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1595 <<
"No scaling method found for " << bc.blockName;
1599 double val =
scale * bc.val;
1602 int nb_integration_pts = OP::getGaussPts().size2();
1603 int nb_base_functions = data.
getN().size2();
1605 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1608 t_N(
i) = t_normal(
i);
1612 t_P(
i,
j) = t_N(
i) * t_N(
j);
1617 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
1621 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
1623 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1625 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1626 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
1630 for (; bb != nb_base_functions; ++bb)
1644template <AssemblyType A>
1650 double time = OP::getFEMethod()->ts_t;
1655 int row_nb_dofs = row_data.
getIndices().size();
1656 int col_nb_dofs = col_data.
getIndices().size();
1657 auto &locMat = OP::locMat;
1658 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1664 for (
auto &bc : (*bcDispPtr)) {
1666 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1668 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1669 auto t_w = OP::getFTensor0IntegrationWeight();
1675 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1676 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1679 <<
"No scaling method found for " << bc.blockName;
1682 int nb_integration_pts = OP::getGaussPts().size2();
1683 int row_nb_dofs = row_data.
getIndices().size();
1684 int col_nb_dofs = col_data.
getIndices().size();
1685 int nb_base_functions = row_data.
getN().size2();
1690 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1693 t_N(
i) = t_normal(
i);
1697 t_P(
i,
j) = t_N(
i) * t_N(
j);
1700 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1703 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1704 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1707 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1708 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1715 for (; rr != nb_base_functions; ++rr)
1722 locMat *= OP::getMeasure();
1729template <AssemblyType A>
1735 double time = OP::getFEMethod()->ts_t;
1740 int row_nb_dofs = row_data.
getIndices().size();
1741 int col_nb_dofs = col_data.
getIndices().size();
1742 auto &locMat = OP::locMat;
1743 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1749 for (
auto &bc : (*bcDispPtr)) {
1751 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1753 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1754 auto t_w = OP::getFTensor0IntegrationWeight();
1763 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1764 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1767 <<
"No scaling method found for " << bc.blockName;
1770 int nb_integration_pts = OP::getGaussPts().size2();
1771 int nb_base_functions = row_data.
getN().size2();
1774 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1777 t_N(
i) = t_normal(
i);
1781 t_P(
i,
j) = t_N(
i) * t_N(
j);
1786 t_d_res(
i,
j) = t_Q(
i,
j);
1789 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1790 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1793 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1795 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1802 for (; rr != nb_base_functions; ++rr)
1809 locMat *= OP::getMeasure();
1816 return OP::iNtegrate(data);
1821 return OP::iNtegrate(row_data, col_data);
1826 return OP::iNtegrate(row_data, col_data);
1829template <AssemblyType A>
1833 double time = OP::getFEMethod()->ts_t;
1841 for (
auto &bc : (*bcDispPtr)) {
1843 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1848 auto [block_name, v_analytical_expr] =
1853 int nb_integration_pts = OP::getGaussPts().size2();
1854 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1855 auto t_w = OP::getFTensor0IntegrationWeight();
1856 int nb_base_functions = data.
getN().size2() / 3;
1858 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1864 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1866 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1867 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1870 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1872 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1876 for (; bb != nb_base_functions; ++bb)
1890 return OP::iNtegrate(data);
1899 int nb_integration_pts = getGaussPts().size2();
1900 int nb_base_functions = data.
getN().size2();
1902 double time = getFEMethod()->ts_t;
1908 if (this->locF.size() != nb_dofs)
1910 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1913 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1916 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1918 auto t_w = getFTensor0IntegrationWeight();
1919 auto t_coords = getFTensor1CoordsAtGaussPts();
1921 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1923 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1925 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1926 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1928 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1929 t_f(
i) -= (time_scale * t_w * t_row_base * tau) * (t_val(
i) *
scale);
1934 for (; rr != nb_base_functions; ++rr)
1939 this->locF *= getMeasure();
1945 for (
auto &bc : *(bcData)) {
1946 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1948 double time_scale = 1;
1949 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1950 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1956 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1960 return -y * (y - 1) / 0.25;
1962 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1965 bc, [](
double,
double,
double) {
return 1; }, time_scale);
1979 int nb_integration_pts = getGaussPts().size2();
1980 int nb_base_functions = data.
getN().size2();
1982 double time = getFEMethod()->ts_t;
1988 if (this->locF.size() != nb_dofs)
1990 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1993 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1998 auto t_w = getFTensor0IntegrationWeight();
1999 auto t_coords = getFTensor1CoordsAtGaussPts();
2000 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2001 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2003 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2005 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2007 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2019 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2022 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2023 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2024 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2026 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2028 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2030 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2032 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2033 t_f(
i) += t_row_base * t_val(
i);
2038 for (; rr != nb_base_functions; ++rr)
2053 for (
auto &bc : *(bcData)) {
2054 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2056 double time_scale = 1;
2057 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2058 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2064 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2071template <AssemblyType A>
2082 double time = OP::getFEMethod()->ts_t;
2087 int nb_base_functions = row_data.
getN().size2();
2088 int row_nb_dofs = row_data.
getIndices().size();
2089 int col_nb_dofs = col_data.
getIndices().size();
2090 int nb_integration_pts = OP::getGaussPts().size2();
2091 auto &locMat = OP::locMat;
2092 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2095auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2100 auto t_w = OP::getFTensor0IntegrationWeight();
2101 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2102 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2103 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2105 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2108 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2118 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2119 auto t_val = time_scale * t_w * tau * val;
2122 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2123 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2126 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2128 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2129 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2130 t_kd(
j,
l) * t_diff_col_base(N0)
2134 (FTensor::levi_civita<double>(
i,
j,
k) *
2135 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2136 t_kd(
k,
l) * t_diff_col_base(N1);
2138 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2145 for (; rr != nb_base_functions; ++rr)
2161 for (
auto &bc : *(bcData)) {
2162 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2164 double time_scale = 1;
2165 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2166 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2172 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2183 return OP::iNtegrate(row_data, col_data);
2193 int nb_integration_pts = getGaussPts().size2();
2194 int nb_base_functions = data.
getN().size2();
2196 double time = getFEMethod()->ts_t;
2202 if (this->locF.size() != nb_dofs)
2204 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2209 for (
auto &bc : *(bcData)) {
2210 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2214 auto [block_name, v_analytical_expr] =
2217 getFTensor1FromMat<3>(v_analytical_expr);
2219 auto t_w = getFTensor0IntegrationWeight();
2220 auto t_coords = getFTensor1CoordsAtGaussPts();
2222 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2224 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2226 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2228 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2229 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2234 for (; rr != nb_base_functions; ++rr)
2240 this->locF *= getMeasure();
2249 int nb_integration_pts = row_data.
getN().size1();
2250 int row_nb_dofs = row_data.
getIndices().size();
2251 int col_nb_dofs = col_data.
getIndices().size();
2252 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2254 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2257 auto v = getVolume();
2258 auto t_w = getFTensor0IntegrationWeight();
2259 int row_nb_base_functions = row_data.
getN().size2();
2261 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2264 for (; rr != row_nb_dofs / 3; ++rr) {
2266 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2267 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2268 double div_col_base = t_col_diff_base_fun(
i,
i);
2269 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2271 ++t_col_diff_base_fun;
2275 for (; rr != row_nb_base_functions; ++rr)
2286 if (alphaW < std::numeric_limits<double>::epsilon() &&
2287 alphaRho < std::numeric_limits<double>::epsilon())
2290 const int nb_integration_pts = row_data.
getN().size1();
2291 const int row_nb_dofs = row_data.
getIndices().size();
2292 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2294 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2300 auto v = getVolume();
2301 auto t_w = getFTensor0IntegrationWeight();
2303 auto piola_scale = dataAtPts->piolaScale;
2304 auto alpha_w = alphaW / piola_scale;
2305 auto alpha_rho = alphaRho / piola_scale;
2307 int row_nb_base_functions = row_data.
getN().size2();
2310 double ts_scale = alpha_w * getTSa();
2311 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2312 ts_scale += alpha_rho * getTSaa();
2314 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2315 double a =
v * t_w * ts_scale;
2318 for (; rr != row_nb_dofs / 3; ++rr) {
2321 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2322 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2323 const double b =
a * t_row_base_fun * t_col_base_fun;
2332 for (; rr != row_nb_base_functions; ++rr)
2351 int nb_integration_pts = row_data.
getN().size1();
2352 int row_nb_dofs = row_data.
getIndices().size();
2353 int col_nb_dofs = col_data.
getIndices().size();
2354 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2357 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2359 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2361 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2363 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2365 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2367 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2370 auto v = getVolume();
2371 auto t_w = getFTensor0IntegrationWeight();
2373 int row_nb_base_functions = row_data.
getN().size2();
2376 auto t_approx_P_adjoint_log_du_dP =
2377 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2379 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2382 for (; rr != row_nb_dofs / 6; ++rr) {
2385 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2387 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2389 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2397 for (; rr != row_nb_base_functions; ++rr)
2400 ++t_approx_P_adjoint_log_du_dP;
2416 int nb_integration_pts = row_data.
getN().size1();
2417 int row_nb_dofs = row_data.
getIndices().size();
2418 int col_nb_dofs = col_data.
getIndices().size();
2419 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2421 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2425 auto v = getVolume();
2426 auto t_w = getFTensor0IntegrationWeight();
2429 int row_nb_base_functions = row_data.
getN().size2();
2431 auto t_approx_P_adjoint_log_du_dP =
2432 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2434 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2437 for (; rr != row_nb_dofs / 6; ++rr) {
2438 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2439 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2440 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2442 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2449 for (; rr != row_nb_base_functions; ++rr)
2452 ++t_approx_P_adjoint_log_du_dP;
2464 int row_nb_dofs = row_data.
getIndices().size();
2465 int col_nb_dofs = col_data.
getIndices().size();
2466 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2469 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2471 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2473 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2475 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2477 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2479 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2489 auto v = getVolume();
2490 auto t_w = getFTensor0IntegrationWeight();
2491 auto t_approx_P_adjoint_log_du_domega =
2492 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2494 int row_nb_base_functions = row_data.
getN().size2();
2497 int nb_integration_pts = row_data.
getN().size1();
2499 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2503 for (; rr != row_nb_dofs / 6; ++rr) {
2505 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2506 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2507 double v =
a * t_row_base_fun * t_col_base_fun;
2508 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2515 for (; rr != row_nb_base_functions; ++rr)
2519 ++t_approx_P_adjoint_log_du_domega;
2528 int nb_integration_pts = getGaussPts().size2();
2529 int row_nb_dofs = row_data.
getIndices().size();
2530 int col_nb_dofs = col_data.
getIndices().size();
2531 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2535 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2536 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2538 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2539 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2541 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2542 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2550 auto v = getVolume();
2551 auto t_w = getFTensor0IntegrationWeight();
2552 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2553 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2554 int row_nb_base_functions = row_data.
getN().size2();
2556 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2559 for (; rr != row_nb_dofs / 3; ++rr) {
2560 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2561 const double b =
a * t_row_base_fun;
2563 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2564 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2570 for (; rr != row_nb_base_functions; ++rr) {
2574 ++t_levi_kirchhoff_du;
2590 int nb_integration_pts = getGaussPts().size2();
2591 int row_nb_dofs = row_data.
getIndices().size();
2592 int col_nb_dofs = col_data.
getIndices().size();
2593 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2596 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2598 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2600 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2605 auto v = getVolume();
2606 auto t_w = getFTensor0IntegrationWeight();
2608 int row_nb_base_functions = row_data.
getN().size2();
2610 auto t_levi_kirchhoff_dP =
2611 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2613 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2616 for (; rr != row_nb_dofs / 3; ++rr) {
2617 double b =
a * t_row_base_fun;
2619 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2620 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2621 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2627 for (; rr != row_nb_base_functions; ++rr) {
2632 ++t_levi_kirchhoff_dP;
2640 int nb_integration_pts = getGaussPts().size2();
2641 int row_nb_dofs = row_data.
getIndices().size();
2642 int col_nb_dofs = col_data.
getIndices().size();
2644 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2646 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2653 auto v = getVolume();
2654 auto t_w = getFTensor0IntegrationWeight();
2655 auto t_levi_kirchoff_dP =
2656 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2658 int row_nb_base_functions = row_data.
getN().size2();
2661 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2664 for (; rr != row_nb_dofs / 3; ++rr) {
2665 double b =
a * t_row_base_fun;
2666 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2667 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2668 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2669 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2676 for (; rr != row_nb_base_functions; ++rr) {
2680 ++t_levi_kirchoff_dP;
2688 int nb_integration_pts = getGaussPts().size2();
2689 int row_nb_dofs = row_data.
getIndices().size();
2690 int col_nb_dofs = col_data.
getIndices().size();
2691 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2694 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2696 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2698 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2711 auto v = getVolume();
2712 auto ts_a = getTSa();
2713 auto t_w = getFTensor0IntegrationWeight();
2714 auto t_levi_kirchhoff_domega =
2715 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2716 int row_nb_base_functions = row_data.
getN().size2();
2721 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2723 double c = (
a * alphaOmega) * (ts_a );
2726 for (; rr != row_nb_dofs / 3; ++rr) {
2727 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2728 const double b =
a * t_row_base_fun;
2731 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2732 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2733 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2741 for (; rr != row_nb_base_functions; ++rr) {
2746 ++t_levi_kirchhoff_domega;
2754 if (dataAtPts->matInvD.size1() ==
size_symm &&
2755 dataAtPts->matInvD.size2() ==
size_symm) {
2768 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2771 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2773 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2775 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2780 int nb_integration_pts = getGaussPts().size2();
2781 int row_nb_dofs = row_data.
getIndices().size();
2782 int col_nb_dofs = col_data.
getIndices().size();
2784 auto v = getVolume();
2785 auto t_w = getFTensor0IntegrationWeight();
2786 int row_nb_base_functions = row_data.
getN().size2() / 3;
2793 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2794 &*dataAtPts->matInvD.data().begin());
2797 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2801 for (; rr != row_nb_dofs / 3; ++rr) {
2803 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2804 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2805 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2813 for (; rr != row_nb_base_functions; ++rr)
2826 if (dataAtPts->matInvD.size1() ==
size_symm &&
2827 dataAtPts->matInvD.size2() ==
size_symm) {
2841 int nb_integration_pts = getGaussPts().size2();
2842 int row_nb_dofs = row_data.
getIndices().size();
2843 int col_nb_dofs = col_data.
getIndices().size();
2845 auto v = getVolume();
2846 auto t_w = getFTensor0IntegrationWeight();
2847 int row_nb_base_functions = row_data.
getN().size2() / 9;
2854 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2855 &*dataAtPts->matInvD.data().begin());
2858 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2862 for (; rr != row_nb_dofs; ++rr) {
2864 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2866 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2873 for (; rr != row_nb_base_functions; ++rr)
2884 if (dataAtPts->matInvD.size1() ==
size_symm &&
2885 dataAtPts->matInvD.size2() ==
size_symm) {
2899 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2902 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
2907 int nb_integration_pts = getGaussPts().size2();
2908 int row_nb_dofs = row_data.
getIndices().size();
2909 int col_nb_dofs = col_data.
getIndices().size();
2911 auto v = getVolume();
2912 auto t_w = getFTensor0IntegrationWeight();
2913 int row_nb_base_functions = row_data.
getN().size2() / 9;
2922 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2923 &*dataAtPts->matInvD.data().begin());
2926 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2929 auto t_m = get_ftensor1(
K, 0, 0);
2932 for (; rr != row_nb_dofs; ++rr) {
2934 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2935 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2943 for (; rr != row_nb_base_functions; ++rr)
2962 int nb_integration_pts = row_data.
getN().size1();
2963 int row_nb_dofs = row_data.
getIndices().size();
2964 int col_nb_dofs = col_data.
getIndices().size();
2966 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2969 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2971 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2973 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2978 auto v = getVolume();
2979 auto t_w = getFTensor0IntegrationWeight();
2980 int row_nb_base_functions = row_data.
getN().size2() / 3;
2983 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2985 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2989 for (; rr != row_nb_dofs / 3; ++rr) {
2992 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2995 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2996 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2997 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
3005 for (; rr != row_nb_base_functions; ++rr)
3025 int nb_integration_pts = row_data.
getN().size1();
3026 int row_nb_dofs = row_data.
getIndices().size();
3027 int col_nb_dofs = col_data.
getIndices().size();
3029 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3031 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3034 auto v = getVolume();
3035 auto t_w = getFTensor0IntegrationWeight();
3036 int row_nb_base_functions = row_data.
getN().size2() / 9;
3039 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3040 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3044 for (; rr != row_nb_dofs; ++rr) {
3047 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3050 auto t_m = get_ftensor2(
K, rr, 0);
3051 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3052 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3060 for (; rr != row_nb_base_functions; ++rr)
3073 if (tagSense != getSkeletonSense())
3076 auto get_tag = [&](
auto name) {
3077 auto &mob = getPtrFE()->mField.get_moab();
3083 auto get_tag_value = [&](
auto &&tag,
int dim) {
3084 auto &mob = getPtrFE()->mField.get_moab();
3085 auto face = getSidePtrFE()->getFEEntityHandle();
3086 std::vector<double> value(dim);
3087 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3091 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3092 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3094 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3095 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3100 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3101 Tag th_detF = create_tag(
"detF", 1);
3102 Tag th_traction = create_tag(
"traction", 3);
3103 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3105 Tag th_energy = create_tag(
"Energy", 1);
3107 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3108 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3109 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3111 auto t_normal = getFTensor1NormalsAtGaussPts();
3112 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3124 if (dataAtPts->energyAtPts.size() == 0) {
3126 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3127 dataAtPts->energyAtPts.clear();
3136 auto set_float_precision = [](
const double x) {
3137 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3144 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3146 v = set_float_precision(
v);
3147 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3154 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3157 for (
auto &
a :
v.data())
3158 a = set_float_precision(
a);
3159 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3160 &*
v.data().begin());
3168 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3170 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3172 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3174 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3176 t_m(
i,
j) = t_d(
i,
j);
3177 for (
auto &
v :
m.data())
3178 v = set_float_precision(
v);
3179 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3180 &*
m.data().begin());
3184 const auto nb_gauss_pts = getGaussPts().size2();
3185 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3188 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3190 t_traction(
i) *= tagSense;
3191 CHKERR save_vec_tag(th_traction, t_traction, gg);
3193 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3194 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3195 CHKERR save_scal_tag(th_energy, t_energy, gg);
3199 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3200 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3201 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3210 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3211 std::vector<FieldSpace> spaces, std::string geom_field_name,
3212 boost::shared_ptr<Range> crack_front_edges_ptr) {
3215 constexpr bool scale_l2 =
false;
3219 "Scale L2 Ainsworth Legendre base is not implemented");
3228 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3229 std::vector<FieldSpace> spaces, std::string geom_field_name,
3230 boost::shared_ptr<Range> crack_front_edges_ptr) {
3233 constexpr bool scale_l2 =
false;
3237 "Scale L2 Ainsworth Legendre base is not implemented");
3246 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3247 std::vector<FieldSpace> spaces, std::string geom_field_name,
3248 boost::shared_ptr<Range> crack_front_edges_ptr,
3249 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3250 boost::shared_ptr<MatrixDouble> inv_jac) {
3254 auto jac = boost::make_shared<MatrixDouble>();
3255 auto det = boost::make_shared<VectorDouble>();
3257 geom_field_name, jac));
3263 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3265 if (scale_l2_ainsworth_legendre_base) {
3273 boost::shared_ptr<MatrixDouble> jac,
3274 boost::shared_ptr<Range> edges_ptr)
3283 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3286 return OP::doWork(side, type, data);
3291 boost::shared_ptr<Range> edgesPtr;
3294 auto jac = boost::make_shared<MatrixDouble>();
3295 auto det = boost::make_shared<VectorDouble>();
3297 geom_field_name, jac,
3369 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3370 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3371 if (getNinTheLoop() == 0) {
3372 dataAtPts->faceMaterialForceAtPts.clear();
3373 dataAtPts->normalPressureAtPts.clear();
3375 auto loop_size = getLoopSize();
3376 if (loop_size == 1) {
3377 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3378 auto pstatus = numebered_fe_ptr->getPStatus();
3379 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3386 auto t_normal = getFTensor1NormalsAtGaussPts();
3387 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3388 dataAtPts->faceMaterialForceAtPts);
3391 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3395 getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
3396 auto t_grad_u_gamma =
3397 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3399 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3400 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
3422 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3423 t_N(
I) = t_normal(
I);
3426 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3428 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3431 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3434 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3436 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3442 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3445 t_N(
I) = t_normal(
I);
3450 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3453 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3454 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3457 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3460 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3463 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3471 "Grffith energy release "
3472 "selector not implemented");
3476 auto side_fe_ptr = getSidePtrFE();
3477 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3478 auto pstatus = side_fe_mi_ptr->getPStatus();
3480 auto owner = side_fe_mi_ptr->getOwnerProc();
3482 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3483 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3484 << getNinTheLoop() <<
" loop size " << getLoopSize();
3496 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3497 auto pstatus = fe_mi_ptr->getPStatus();
3499 auto owner = fe_mi_ptr->getOwnerProc();
3501 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3502 <<
" " << getPtrFE()->mField.get_comm_rank();
3510 double face_pressure = 0.;
3511 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3512 dataAtPts->faceMaterialForceAtPts);
3515 auto t_w = getFTensor0IntegrationWeight();
3516 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3517 t_face_T(
I) += t_w * t_T(
I);
3518 face_pressure += t_w * t_p;
3523 t_face_T(
I) *= getMeasure();
3524 face_pressure *= getMeasure();
3526 auto get_tag = [&](
auto name,
auto dim) {
3527 auto &moab = getPtrFE()->mField.get_moab();
3529 double def_val[] = {0., 0., 0.};
3530 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3531 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3536 auto set_tag = [&](
auto &&tag,
auto ptr) {
3537 auto &moab = getPtrFE()->mField.get_moab();
3538 auto face = getPtrFE()->getFEEntityHandle();
3539 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3542 set_tag(
get_tag(
"MaterialForce", 3), &t_face_T(0));
3543 set_tag(
get_tag(
"FacePressure", 1), &face_pressure);
3549template <
typename OP_PTR>
3551 const std::string block_name) {
3553 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
3555 auto ts_time = op_ptr->getTStime();
3556 auto ts_time_step = op_ptr->getTStimeStep();
3562 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
3563 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
3565 auto v_analytical_expr =
3567 m_ref_coords, m_ref_normals, block_name);
3570 if (v_analytical_expr.size2() != nb_gauss_pts)
3572 "Wrong number of integration pts");
3575 return std::make_tuple(block_name, v_analytical_expr);
3579#ifdef ENABLE_PYTHON_BINDING
3581 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3586 auto main_module = bp::import(
"__main__");
3587 mainNamespace = main_module.attr(
"__dict__");
3588 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3590 analyticalDispFun = mainNamespace[
"analytical_disp"];
3591 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3592 analyticalExternalStrainFun = mainNamespace[
"analytical_external_strain"];
3594 }
catch (bp::error_already_set
const &) {
3604 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3605 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3606 const std::string &block_name, np::ndarray &analytical_expr
3613 analytical_expr = bp::extract<np::ndarray>(
3614 analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3616 }
catch (bp::error_already_set
const &) {
3626 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3627 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3628 const std::string& block_name, np::ndarray &analytical_expr
3635 analytical_expr = bp::extract<np::ndarray>(
3636 analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3638 }
catch (bp::error_already_set
const &) {
3646 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalExternalStrain(
3648 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3649 const std::string& block_name, np::ndarray &analytical_expr
3656 analytical_expr = bp::extract<np::ndarray>(
3657 analyticalExternalStrainFun(delta_t,
t, x, y, z, block_name));
3659 }
catch (bp::error_already_set
const &) {
3667boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3669inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3671 auto dtype = np::dtype::get_builtin<double>();
3672 auto size = bp::make_tuple(nb_gauss_pts);
3673 auto stride = bp::make_tuple(3 *
sizeof(
double));
3674 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3683 const std::string block_name) {
3685#ifdef ENABLE_PYTHON_BINDING
3686 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3691 bp::list python_coords;
3692 bp::list python_normals;
3694 for (
int idx = 0; idx < 3; ++idx) {
3695 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3696 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3699 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3700 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3702 std::regex reg_disp_name(disp_block_name);
3703 std::regex reg_traction_name(traction_block_name);
3705 np::ndarray np_analytical_expr = np::empty(
3706 bp::make_tuple(nb_gauss_pts, 3), np::dtype::get_builtin<double>());
3708 if (std::regex_match(block_name, reg_disp_name)) {
3710 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3711 bp::extract<np::ndarray>(python_coords[1]),
3712 bp::extract<np::ndarray>(python_coords[2]),
3713 bp::extract<np::ndarray>(python_normals[0]),
3714 bp::extract<np::ndarray>(python_normals[1]),
3715 bp::extract<np::ndarray>(python_normals[2]),
3716 block_name, np_analytical_expr),
3717 "Failed analytical_disp() python call");
3718 }
else if (std::regex_match(block_name, reg_traction_name)) {
3720 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3721 bp::extract<np::ndarray>(python_coords[1]),
3722 bp::extract<np::ndarray>(python_coords[2]),
3723 bp::extract<np::ndarray>(python_normals[0]),
3724 bp::extract<np::ndarray>(python_normals[1]),
3725 bp::extract<np::ndarray>(python_normals[2]),
3726 block_name, np_analytical_expr),
3727 "Failed analytical_traction() python call");
3730 "Unknown analytical block");
3733 double *analytical_expr_val_ptr =
3734 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3737 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3738 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3739 for (
int idx = 0; idx < 3; ++idx)
3740 v_analytical_expr(idx, gg) =
3741 *(analytical_expr_val_ptr + (3 * gg + idx));
3743 return v_analytical_expr;
3754 const std::string block_name) {
3756#ifdef ENABLE_PYTHON_BINDING
3757 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3761 bp::list python_coords;
3763 for (
int idx = 0; idx < 3; ++idx) {
3764 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3767 auto externalstrain_block_name =
"(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
3769 std::regex reg_externalstrain_name(externalstrain_block_name);
3771 np::ndarray np_analytical_expr = np::empty(
3772 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
3774 if (std::regex_match(block_name, reg_externalstrain_name)) {
3775 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalExternalStrain(
3776 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3777 bp::extract<np::ndarray>(python_coords[1]),
3778 bp::extract<np::ndarray>(python_coords[2]), block_name,
3779 np_analytical_expr),
3780 "Failed analytical_external_strain() python call");
3783 "Unknown analytical block");
3786 double *analytical_expr_val_ptr =
3787 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3790 v_analytical_expr.resize(nb_gauss_pts,
false);
3791 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg)
3792 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
3794 return v_analytical_expr;
3803 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3804 boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
3805 boost::shared_ptr<Range> ents_ptr)
3806 :
OP(broken_base_side_data, ents_ptr) {
3807 this->sourceVec = vec;
3808 this->betaCoeff = beta_coeff;
3812 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3813 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3814 :
OP(broken_base_side_data, ents_ptr) {
3815 this->sourceVec = boost::shared_ptr<MatrixDouble>();
3816 this->betaCoeff = beta_coeff;
3825 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3830 if (!brokenBaseSideData) {
3835 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
3843 OP::nbIntegrationPts = OP::getGaussPts().size2();
3845 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
3847 OP::locF.resize(OP::nbRows,
false);
3852 CHKERR this->aSsemble(row_data);
3856 switch (OP::opType) {
3858 for (
auto &bd : *brokenBaseSideData) {
3860 boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
3861 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
3862 this->sourceVec.reset();
3867 (std::string(
"wrong op type ") +
3868 OpBaseDerivativesBase::OpTypeNames[OP::opType])
3876 const std::string row_field,
3877 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3878 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3880 brokenBaseSideDataPtr(broken_base_side_data) {
3881 this->betaCoeff = beta_coeff;
3887 for (
auto &bd : (*brokenBaseSideDataPtr)) {
3892 if (this->sourceVec->size1() !=
SPACE_DIM) {
3894 "Inconsistent size of the source vector");
3896 if (this->sourceVec->size2() != OP::getGaussPts().size2()) {
3898 "Inconsistent size of the source vector");
3902 CHKERR OP::iNtegrate(data);
3904 this->sourceVec.reset();
3910 std::string row_field,
3911 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3912 ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
3913 boost::shared_ptr<Range> ents_ptr)
3914 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
3916 this->betaCoeff = beta;
3921 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3922 ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
3923 :
OP(broken_base_side_data, ents_ptr) {
3925 this->betaCoeff = beta;
3934 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3939 if (!brokenBaseSideData) {
3944 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
3945 EntityType col_type,
3950 auto check_if_assemble_transpose = [&] {
3952 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
3956 }
else if (OP::assembleTranspose) {
3962 OP::rowSide = row_side;
3963 OP::rowType = row_type;
3964 OP::colSide = col_side;
3965 OP::colType = col_type;
3967 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
3970 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
3974 switch (OP::opType) {
3977 for (
auto &bd : *brokenBaseSideData) {
3980 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
3982 "base functions not set");
3986 OP::nbRows = bd.getData().getIndices().size();
3989 OP::nbIntegrationPts = OP::getGaussPts().size2();
3990 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
3998 bd.getSide(), bd.getSide(),
4001 bd.getType(), bd.getType(),
4004 bd.getData(), bd.getData()
4013 (std::string(
"wrong op type ") +
4014 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
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
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)
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)
