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);
1317 : getFEMethod()->ts_t;
1320 double scale = scalingMethodPtr->getScale(time);
1325 int nb_base_functions = data.
getN().size2();
1327 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1329 auto t_nf = get_ftensor2(nF);
1332 t_symm_stress(
i,
j) =
1333 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1336 t_residual(
L) = t_L(
i,
j,
L) * (
scale * t_symm_stress(
i,
j));
1339 for (; bb != nb_dofs / 6; ++bb) {
1340 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1344 for (; bb != nb_base_functions; ++bb)
1348 ++t_internal_stress;
1358 int nb_integration_pts = data.
getN().size1();
1359 auto v = getVolume();
1360 auto t_w = getFTensor0IntegrationWeight();
1362 auto get_ftensor2 = [](
auto &
v) {
1364 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1367 auto t_internal_stress =
1368 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1377 : getFEMethod()->ts_t;
1380 double scale = scalingMethodPtr->getScale(time);
1382 int nb_base_functions = data.
getN().size2();
1384 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1386 auto t_nf = get_ftensor2(nF);
1389 t_residual(
L) = t_L(M,
L) * (
scale * t_internal_stress(M));
1392 for (; bb != nb_dofs / 6; ++bb) {
1393 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1397 for (; bb != nb_base_functions; ++bb)
1401 ++t_internal_stress;
1406template <AssemblyType A>
1412 for (
auto &bc : (*bcDispPtr)) {
1414 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1417 int nb_integration_pts = OP::getGaussPts().size2();
1418 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1419 auto t_w = OP::getFTensor0IntegrationWeight();
1420 int nb_base_functions = data.
getN().size2() / 3;
1427 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1429 scale *= scalingMethodsMap.at(bc.blockName)
1432 scale *= scalingMethodsMap.at(bc.blockName)
1433 ->getScale(OP::getFEMethod()->ts_t);
1437 <<
"No scaling method found for " << bc.blockName;
1444 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1445 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1447 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1449 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1453 for (; bb != nb_base_functions; ++bb)
1465 return OP::iNtegrate(data);
1468template <AssemblyType A>
1476 double time = OP::getFEMethod()->ts_t;
1484 for (
auto &bc : (*bcRotPtr)) {
1486 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1488 int nb_integration_pts = OP::getGaussPts().size2();
1489 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1490 auto t_w = OP::getFTensor0IntegrationWeight();
1492 int nb_base_functions = data.
getN().size2() / 3;
1495 auto get_ftensor1 = [](
auto &
v) {
1508 auto get_rotation_angle = [&]() {
1509 double theta = bc.theta;
1510 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1511 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1516 auto get_rotation = [&](
auto theta) {
1518 if (bc.vals.size() == 7) {
1519 t_omega(0) = bc.vals[4];
1520 t_omega(1) = bc.vals[5];
1521 t_omega(2) = bc.vals[6];
1524 t_omega(
i) = OP::getFTensor1Normal()(
i);
1527 t_omega(
i) *= theta;
1534 auto t_R = get_rotation(get_rotation_angle());
1535 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1537 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1539 t_delta(
i) = t_center(
i) - t_coords(
i);
1541 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1543 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1545 for (; bb != nb_dofs / 3; ++bb) {
1546 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1550 for (; bb != nb_base_functions; ++bb)
1563 return OP::iNtegrate(data);
1566template <AssemblyType A>
1570 double time = OP::getFEMethod()->ts_t;
1578 for (
auto &bc : (*bcDispPtr)) {
1580 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1582 for (
auto &bd : (*brokenBaseSideDataPtr)) {
1584 auto t_approx_P = getFTensor2FromMat<3, 3>(bd.getFlux());
1585 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1586 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1587 auto t_w = OP::getFTensor0IntegrationWeight();
1595 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1596 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1599 <<
"No scaling method found for " << bc.blockName;
1603 double val =
scale * bc.val;
1606 int nb_integration_pts = OP::getGaussPts().size2();
1607 int nb_base_functions = data.
getN().size2();
1609 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1612 t_N(
i) = t_normal(
i);
1616 t_P(
i,
j) = t_N(
i) * t_N(
j);
1621 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
1625 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
1627 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1629 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1630 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
1634 for (; bb != nb_base_functions; ++bb)
1648template <AssemblyType A>
1654 double time = OP::getFEMethod()->ts_t;
1659 int row_nb_dofs = row_data.
getIndices().size();
1660 int col_nb_dofs = col_data.
getIndices().size();
1661 auto &locMat = OP::locMat;
1662 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1668 for (
auto &bc : (*bcDispPtr)) {
1670 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1672 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1673 auto t_w = OP::getFTensor0IntegrationWeight();
1679 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1680 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1683 <<
"No scaling method found for " << bc.blockName;
1686 int nb_integration_pts = OP::getGaussPts().size2();
1687 int row_nb_dofs = row_data.
getIndices().size();
1688 int col_nb_dofs = col_data.
getIndices().size();
1689 int nb_base_functions = row_data.
getN().size2();
1694 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1697 t_N(
i) = t_normal(
i);
1701 t_P(
i,
j) = t_N(
i) * t_N(
j);
1704 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1707 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1708 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1711 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1712 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1719 for (; rr != nb_base_functions; ++rr)
1726 locMat *= OP::getMeasure();
1733template <AssemblyType A>
1739 double time = OP::getFEMethod()->ts_t;
1744 int row_nb_dofs = row_data.
getIndices().size();
1745 int col_nb_dofs = col_data.
getIndices().size();
1746 auto &locMat = OP::locMat;
1747 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1753 for (
auto &bc : (*bcDispPtr)) {
1755 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1757 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1758 auto t_w = OP::getFTensor0IntegrationWeight();
1767 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1768 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1771 <<
"No scaling method found for " << bc.blockName;
1774 int nb_integration_pts = OP::getGaussPts().size2();
1775 int nb_base_functions = row_data.
getN().size2();
1778 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1781 t_N(
i) = t_normal(
i);
1785 t_P(
i,
j) = t_N(
i) * t_N(
j);
1790 t_d_res(
i,
j) = t_Q(
i,
j);
1793 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1794 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1797 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1799 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1806 for (; rr != nb_base_functions; ++rr)
1813 locMat *= OP::getMeasure();
1820 return OP::iNtegrate(data);
1825 return OP::iNtegrate(row_data, col_data);
1830 return OP::iNtegrate(row_data, col_data);
1833template <AssemblyType A>
1837 double time = OP::getFEMethod()->ts_t;
1845 for (
auto &bc : (*bcDispPtr)) {
1847 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1852 auto [block_name, v_analytical_expr] =
1857 int nb_integration_pts = OP::getGaussPts().size2();
1858 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1859 auto t_w = OP::getFTensor0IntegrationWeight();
1860 int nb_base_functions = data.
getN().size2() / 3;
1862 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1868 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1870 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1871 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1874 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1876 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1880 for (; bb != nb_base_functions; ++bb)
1894 return OP::iNtegrate(data);
1903 int nb_integration_pts = getGaussPts().size2();
1904 int nb_base_functions = data.
getN().size2();
1906 double time = getFEMethod()->ts_t;
1912 if (this->locF.size() != nb_dofs)
1914 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1917 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1920 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1922 auto t_w = getFTensor0IntegrationWeight();
1923 auto t_coords = getFTensor1CoordsAtGaussPts();
1925 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1927 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1929 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1930 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1932 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1933 t_f(
i) -= (time_scale * t_w * t_row_base * tau) * (t_val(
i) *
scale);
1938 for (; rr != nb_base_functions; ++rr)
1943 this->locF *= getMeasure();
1949 for (
auto &bc : *(bcData)) {
1950 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1952 double time_scale = 1;
1953 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1954 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1960 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1964 return -y * (y - 1) / 0.25;
1966 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1969 bc, [](
double,
double,
double) {
return 1; }, time_scale);
1983 int nb_integration_pts = getGaussPts().size2();
1984 int nb_base_functions = data.
getN().size2();
1986 double time = getFEMethod()->ts_t;
1992 if (this->locF.size() != nb_dofs)
1994 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1997 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2002 auto t_w = getFTensor0IntegrationWeight();
2003 auto t_coords = getFTensor1CoordsAtGaussPts();
2004 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2005 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2007 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2009 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2011 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2023 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2026 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2027 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2028 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2030 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2032 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2034 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2036 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2037 t_f(
i) += t_row_base * t_val(
i);
2042 for (; rr != nb_base_functions; ++rr)
2057 for (
auto &bc : *(bcData)) {
2058 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2060 double time_scale = 1;
2061 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2062 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2068 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2075template <AssemblyType A>
2086 double time = OP::getFEMethod()->ts_t;
2091 int nb_base_functions = row_data.
getN().size2();
2092 int row_nb_dofs = row_data.
getIndices().size();
2093 int col_nb_dofs = col_data.
getIndices().size();
2094 int nb_integration_pts = OP::getGaussPts().size2();
2095 auto &locMat = OP::locMat;
2096 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2099auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2104 auto t_w = OP::getFTensor0IntegrationWeight();
2105 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2106 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2107 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2109 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2112 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2122 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2123 auto t_val = time_scale * t_w * tau * val;
2126 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2127 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2130 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2132 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2133 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2134 t_kd(
j,
l) * t_diff_col_base(N0)
2138 (FTensor::levi_civita<double>(
i,
j,
k) *
2139 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2140 t_kd(
k,
l) * t_diff_col_base(N1);
2142 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2149 for (; rr != nb_base_functions; ++rr)
2165 for (
auto &bc : *(bcData)) {
2166 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2168 double time_scale = 1;
2169 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2170 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2176 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2187 return OP::iNtegrate(row_data, col_data);
2197 int nb_integration_pts = getGaussPts().size2();
2198 int nb_base_functions = data.
getN().size2();
2200 double time = getFEMethod()->ts_t;
2206 if (this->locF.size() != nb_dofs)
2208 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2213 for (
auto &bc : *(bcData)) {
2214 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2218 auto [block_name, v_analytical_expr] =
2221 getFTensor1FromMat<3>(v_analytical_expr);
2223 auto t_w = getFTensor0IntegrationWeight();
2224 auto t_coords = getFTensor1CoordsAtGaussPts();
2226 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2228 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2230 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2232 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2233 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2238 for (; rr != nb_base_functions; ++rr)
2244 this->locF *= getMeasure();
2253 int nb_integration_pts = row_data.
getN().size1();
2254 int row_nb_dofs = row_data.
getIndices().size();
2255 int col_nb_dofs = col_data.
getIndices().size();
2256 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2258 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2261 auto v = getVolume();
2262 auto t_w = getFTensor0IntegrationWeight();
2263 int row_nb_base_functions = row_data.
getN().size2();
2265 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2268 for (; rr != row_nb_dofs / 3; ++rr) {
2270 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2271 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2272 double div_col_base = t_col_diff_base_fun(
i,
i);
2273 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2275 ++t_col_diff_base_fun;
2279 for (; rr != row_nb_base_functions; ++rr)
2290 if (alphaW < std::numeric_limits<double>::epsilon() &&
2291 alphaRho < std::numeric_limits<double>::epsilon())
2294 const int nb_integration_pts = row_data.
getN().size1();
2295 const int row_nb_dofs = row_data.
getIndices().size();
2296 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2298 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2304 auto v = getVolume();
2305 auto t_w = getFTensor0IntegrationWeight();
2307 auto piola_scale = dataAtPts->piolaScale;
2308 auto alpha_w = alphaW / piola_scale;
2309 auto alpha_rho = alphaRho / piola_scale;
2311 int row_nb_base_functions = row_data.
getN().size2();
2314 double ts_scale = alpha_w * getTSa();
2315 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2316 ts_scale += alpha_rho * getTSaa();
2318 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2319 double a =
v * t_w * ts_scale;
2322 for (; rr != row_nb_dofs / 3; ++rr) {
2325 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2326 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2327 const double b =
a * t_row_base_fun * t_col_base_fun;
2336 for (; rr != row_nb_base_functions; ++rr)
2355 int nb_integration_pts = row_data.
getN().size1();
2356 int row_nb_dofs = row_data.
getIndices().size();
2357 int col_nb_dofs = col_data.
getIndices().size();
2358 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2361 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2363 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2365 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2367 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2369 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2371 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2374 auto v = getVolume();
2375 auto t_w = getFTensor0IntegrationWeight();
2377 int row_nb_base_functions = row_data.
getN().size2();
2380 auto t_approx_P_adjoint_log_du_dP =
2381 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2383 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2386 for (; rr != row_nb_dofs / 6; ++rr) {
2389 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2391 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2393 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2401 for (; rr != row_nb_base_functions; ++rr)
2404 ++t_approx_P_adjoint_log_du_dP;
2420 int nb_integration_pts = row_data.
getN().size1();
2421 int row_nb_dofs = row_data.
getIndices().size();
2422 int col_nb_dofs = col_data.
getIndices().size();
2423 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2425 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2429 auto v = getVolume();
2430 auto t_w = getFTensor0IntegrationWeight();
2433 int row_nb_base_functions = row_data.
getN().size2();
2435 auto t_approx_P_adjoint_log_du_dP =
2436 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2438 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2441 for (; rr != row_nb_dofs / 6; ++rr) {
2442 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2443 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2444 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2446 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2453 for (; rr != row_nb_base_functions; ++rr)
2456 ++t_approx_P_adjoint_log_du_dP;
2468 int row_nb_dofs = row_data.
getIndices().size();
2469 int col_nb_dofs = col_data.
getIndices().size();
2470 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2473 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2475 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2477 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2479 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2481 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2483 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2493 auto v = getVolume();
2494 auto t_w = getFTensor0IntegrationWeight();
2495 auto t_approx_P_adjoint_log_du_domega =
2496 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2498 int row_nb_base_functions = row_data.
getN().size2();
2501 int nb_integration_pts = row_data.
getN().size1();
2503 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2507 for (; rr != row_nb_dofs / 6; ++rr) {
2509 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2510 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2511 double v =
a * t_row_base_fun * t_col_base_fun;
2512 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2519 for (; rr != row_nb_base_functions; ++rr)
2523 ++t_approx_P_adjoint_log_du_domega;
2532 int nb_integration_pts = getGaussPts().size2();
2533 int row_nb_dofs = row_data.
getIndices().size();
2534 int col_nb_dofs = col_data.
getIndices().size();
2535 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2539 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2540 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2542 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2543 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2545 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2546 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2554 auto v = getVolume();
2555 auto t_w = getFTensor0IntegrationWeight();
2556 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2557 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2558 int row_nb_base_functions = row_data.
getN().size2();
2560 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2563 for (; rr != row_nb_dofs / 3; ++rr) {
2564 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2565 const double b =
a * t_row_base_fun;
2567 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2568 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2574 for (; rr != row_nb_base_functions; ++rr) {
2578 ++t_levi_kirchhoff_du;
2594 int nb_integration_pts = getGaussPts().size2();
2595 int row_nb_dofs = row_data.
getIndices().size();
2596 int col_nb_dofs = col_data.
getIndices().size();
2597 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2600 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2602 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2604 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2609 auto v = getVolume();
2610 auto t_w = getFTensor0IntegrationWeight();
2612 int row_nb_base_functions = row_data.
getN().size2();
2614 auto t_levi_kirchhoff_dP =
2615 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2617 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2620 for (; rr != row_nb_dofs / 3; ++rr) {
2621 double b =
a * t_row_base_fun;
2623 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2624 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2625 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2631 for (; rr != row_nb_base_functions; ++rr) {
2636 ++t_levi_kirchhoff_dP;
2644 int nb_integration_pts = getGaussPts().size2();
2645 int row_nb_dofs = row_data.
getIndices().size();
2646 int col_nb_dofs = col_data.
getIndices().size();
2648 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2650 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2657 auto v = getVolume();
2658 auto t_w = getFTensor0IntegrationWeight();
2659 auto t_levi_kirchoff_dP =
2660 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2662 int row_nb_base_functions = row_data.
getN().size2();
2665 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2668 for (; rr != row_nb_dofs / 3; ++rr) {
2669 double b =
a * t_row_base_fun;
2670 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2671 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2672 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2673 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2680 for (; rr != row_nb_base_functions; ++rr) {
2684 ++t_levi_kirchoff_dP;
2692 int nb_integration_pts = getGaussPts().size2();
2693 int row_nb_dofs = row_data.
getIndices().size();
2694 int col_nb_dofs = col_data.
getIndices().size();
2695 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2698 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2700 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2702 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2715 auto v = getVolume();
2716 auto ts_a = getTSa();
2717 auto t_w = getFTensor0IntegrationWeight();
2718 auto t_levi_kirchhoff_domega =
2719 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2720 int row_nb_base_functions = row_data.
getN().size2();
2725 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2727 double c = (
a * alphaOmega) * (ts_a );
2730 for (; rr != row_nb_dofs / 3; ++rr) {
2731 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2732 const double b =
a * t_row_base_fun;
2735 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2736 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2737 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2745 for (; rr != row_nb_base_functions; ++rr) {
2750 ++t_levi_kirchhoff_domega;
2758 if (dataAtPts->matInvD.size1() ==
size_symm &&
2759 dataAtPts->matInvD.size2() ==
size_symm) {
2772 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2775 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2777 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2779 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2784 int nb_integration_pts = getGaussPts().size2();
2785 int row_nb_dofs = row_data.
getIndices().size();
2786 int col_nb_dofs = col_data.
getIndices().size();
2788 auto v = getVolume();
2789 auto t_w = getFTensor0IntegrationWeight();
2790 int row_nb_base_functions = row_data.
getN().size2() / 3;
2797 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2798 &*dataAtPts->matInvD.data().begin());
2801 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2805 for (; rr != row_nb_dofs / 3; ++rr) {
2807 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2808 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2809 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2817 for (; rr != row_nb_base_functions; ++rr)
2830 if (dataAtPts->matInvD.size1() ==
size_symm &&
2831 dataAtPts->matInvD.size2() ==
size_symm) {
2845 int nb_integration_pts = getGaussPts().size2();
2846 int row_nb_dofs = row_data.
getIndices().size();
2847 int col_nb_dofs = col_data.
getIndices().size();
2849 auto v = getVolume();
2850 auto t_w = getFTensor0IntegrationWeight();
2851 int row_nb_base_functions = row_data.
getN().size2() / 9;
2858 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2859 &*dataAtPts->matInvD.data().begin());
2862 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2866 for (; rr != row_nb_dofs; ++rr) {
2868 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2870 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2877 for (; rr != row_nb_base_functions; ++rr)
2888 if (dataAtPts->matInvD.size1() ==
size_symm &&
2889 dataAtPts->matInvD.size2() ==
size_symm) {
2903 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2906 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
2911 int nb_integration_pts = getGaussPts().size2();
2912 int row_nb_dofs = row_data.
getIndices().size();
2913 int col_nb_dofs = col_data.
getIndices().size();
2915 auto v = getVolume();
2916 auto t_w = getFTensor0IntegrationWeight();
2917 int row_nb_base_functions = row_data.
getN().size2() / 9;
2926 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2927 &*dataAtPts->matInvD.data().begin());
2930 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2933 auto t_m = get_ftensor1(
K, 0, 0);
2936 for (; rr != row_nb_dofs; ++rr) {
2938 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2939 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2947 for (; rr != row_nb_base_functions; ++rr)
2966 int nb_integration_pts = row_data.
getN().size1();
2967 int row_nb_dofs = row_data.
getIndices().size();
2968 int col_nb_dofs = col_data.
getIndices().size();
2970 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2973 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2975 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2977 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2982 auto v = getVolume();
2983 auto t_w = getFTensor0IntegrationWeight();
2984 int row_nb_base_functions = row_data.
getN().size2() / 3;
2987 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2989 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2993 for (; rr != row_nb_dofs / 3; ++rr) {
2996 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2999 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3000 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3001 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
3009 for (; rr != row_nb_base_functions; ++rr)
3029 int nb_integration_pts = row_data.
getN().size1();
3030 int row_nb_dofs = row_data.
getIndices().size();
3031 int col_nb_dofs = col_data.
getIndices().size();
3033 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3035 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3038 auto v = getVolume();
3039 auto t_w = getFTensor0IntegrationWeight();
3040 int row_nb_base_functions = row_data.
getN().size2() / 9;
3043 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3044 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3048 for (; rr != row_nb_dofs; ++rr) {
3051 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3054 auto t_m = get_ftensor2(
K, rr, 0);
3055 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3056 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3064 for (; rr != row_nb_base_functions; ++rr)
3077 if (tagSense != getSkeletonSense())
3080 auto get_tag = [&](
auto name) {
3081 auto &mob = getPtrFE()->mField.get_moab();
3087 auto get_tag_value = [&](
auto &&tag,
int dim) {
3088 auto &mob = getPtrFE()->mField.get_moab();
3089 auto face = getSidePtrFE()->getFEEntityHandle();
3090 std::vector<double> value(dim);
3091 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3095 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3096 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3098 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3099 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3104 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3105 Tag th_detF = create_tag(
"detF", 1);
3106 Tag th_traction = create_tag(
"traction", 3);
3107 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3109 Tag th_energy = create_tag(
"Energy", 1);
3111 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3112 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3113 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3115 auto t_normal = getFTensor1NormalsAtGaussPts();
3116 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3128 if (dataAtPts->energyAtPts.size() == 0) {
3130 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3131 dataAtPts->energyAtPts.clear();
3140 auto set_float_precision = [](
const double x) {
3141 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3148 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3150 v = set_float_precision(
v);
3151 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3158 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3161 for (
auto &
a :
v.data())
3162 a = set_float_precision(
a);
3163 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3164 &*
v.data().begin());
3172 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3174 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3176 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3178 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3180 t_m(
i,
j) = t_d(
i,
j);
3181 for (
auto &
v :
m.data())
3182 v = set_float_precision(
v);
3183 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3184 &*
m.data().begin());
3188 const auto nb_gauss_pts = getGaussPts().size2();
3189 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3192 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3194 t_traction(
i) *= tagSense;
3195 CHKERR save_vec_tag(th_traction, t_traction, gg);
3197 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3198 if (!std::isfinite(u_error))
3200 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3201 CHKERR save_scal_tag(th_energy, t_energy, gg);
3205 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3206 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3207 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3216 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3217 std::vector<FieldSpace> spaces, std::string geom_field_name,
3218 boost::shared_ptr<Range> crack_front_edges_ptr) {
3221 constexpr bool scale_l2 =
false;
3225 "Scale L2 Ainsworth Legendre base is not implemented");
3234 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3235 std::vector<FieldSpace> spaces, std::string geom_field_name,
3236 boost::shared_ptr<Range> crack_front_edges_ptr) {
3239 constexpr bool scale_l2 =
false;
3243 "Scale L2 Ainsworth Legendre base is not implemented");
3252 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3253 std::vector<FieldSpace> spaces, std::string geom_field_name,
3254 boost::shared_ptr<Range> crack_front_edges_ptr,
3255 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3256 boost::shared_ptr<MatrixDouble> inv_jac) {
3260 auto jac = boost::make_shared<MatrixDouble>();
3261 auto det = boost::make_shared<VectorDouble>();
3263 geom_field_name, jac));
3269 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3271 if (scale_l2_ainsworth_legendre_base) {
3279 boost::shared_ptr<MatrixDouble> jac,
3280 boost::shared_ptr<Range> edges_ptr)
3289 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3292 return OP::doWork(side, type, data);
3297 boost::shared_ptr<Range> edgesPtr;
3300 auto jac = boost::make_shared<MatrixDouble>();
3301 auto det = boost::make_shared<VectorDouble>();
3303 geom_field_name, jac,
3375 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3376 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3377 if (getNinTheLoop() == 0) {
3378 dataAtPts->faceMaterialForceAtPts.clear();
3379 dataAtPts->normalPressureAtPts.clear();
3381 auto loop_size = getLoopSize();
3382 if (loop_size == 1) {
3383 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3384 auto pstatus = numebered_fe_ptr->getPStatus();
3385 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3392 auto t_normal = getFTensor1NormalsAtGaussPts();
3393 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3394 dataAtPts->faceMaterialForceAtPts);
3397 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3401 getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
3402 auto t_grad_u_gamma =
3403 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3405 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3406 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
3428 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3429 t_N(
I) = t_normal(
I);
3432 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3434 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3437 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3440 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3443 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3450 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3453 t_N(
I) = t_normal(
I);
3458 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3461 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3462 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3465 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3468 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3472 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3481 "Grffith energy release "
3482 "selector not implemented");
3486 auto side_fe_ptr = getSidePtrFE();
3487 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3488 auto pstatus = side_fe_mi_ptr->getPStatus();
3490 auto owner = side_fe_mi_ptr->getOwnerProc();
3492 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3493 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3494 << getNinTheLoop() <<
" loop size " << getLoopSize();
3506 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3507 auto pstatus = fe_mi_ptr->getPStatus();
3509 auto owner = fe_mi_ptr->getOwnerProc();
3511 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3512 <<
" " << getPtrFE()->mField.get_comm_rank();
3520 double face_pressure = 0.;
3521 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3522 dataAtPts->faceMaterialForceAtPts);
3525 auto t_w = getFTensor0IntegrationWeight();
3526 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3527 t_face_T(
I) += t_w * t_T(
I);
3528 face_pressure += t_w * t_p;
3533 t_face_T(
I) *= getMeasure();
3534 face_pressure *= getMeasure();
3536 auto get_tag = [&](
auto name,
auto dim) {
3537 auto &moab = getPtrFE()->mField.get_moab();
3539 double def_val[] = {0., 0., 0.};
3540 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3541 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3546 auto set_tag = [&](
auto &&tag,
auto ptr) {
3547 auto &moab = getPtrFE()->mField.get_moab();
3548 auto face = getPtrFE()->getFEEntityHandle();
3549 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3552 set_tag(
get_tag(
"MaterialForce", 3), &t_face_T(0));
3553 set_tag(
get_tag(
"FacePressure", 1), &face_pressure);
3559template <
typename OP_PTR>
3561 const std::string block_name) {
3563 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
3565 auto ts_time = op_ptr->getTStime();
3566 auto ts_time_step = op_ptr->getTStimeStep();
3572 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
3573 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
3575 auto v_analytical_expr =
3577 m_ref_coords, m_ref_normals, block_name);
3579 if (PetscUnlikely(!v_analytical_expr.size2())) {
3581 "Analytical expression is empty or does not exist, "
3582 "check python file");
3585 return std::make_tuple(block_name, v_analytical_expr);
3589#ifdef ENABLE_PYTHON_BINDING
3591 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3596 auto main_module = bp::import(
"__main__");
3597 mainNamespace = main_module.attr(
"__dict__");
3598 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3600 analyticalDispFun = mainNamespace[
"analytical_disp"];
3601 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3602 analyticalExternalStrainFun = mainNamespace[
"analytical_external_strain"];
3604 }
catch (bp::error_already_set
const &) {
3614 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3615 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3616 const std::string &block_name, np::ndarray &analytical_expr
3623 analytical_expr = bp::extract<np::ndarray>(
3624 analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3626 }
catch (bp::error_already_set
const &) {
3636 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3637 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3638 const std::string& block_name, np::ndarray &analytical_expr
3645 analytical_expr = bp::extract<np::ndarray>(
3646 analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3648 }
catch (bp::error_already_set
const &) {
3656 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalExternalStrain(
3658 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3659 const std::string& block_name, np::ndarray &analytical_expr
3666 analytical_expr = bp::extract<np::ndarray>(
3667 analyticalExternalStrainFun(delta_t,
t, x, y, z, block_name));
3669 }
catch (bp::error_already_set
const &) {
3677boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3679inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3681 auto dtype = np::dtype::get_builtin<double>();
3682 auto size = bp::make_tuple(nb_gauss_pts);
3683 auto stride = bp::make_tuple(3 *
sizeof(
double));
3684 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3693 const std::string block_name) {
3695#ifdef ENABLE_PYTHON_BINDING
3696 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3701 bp::list python_coords;
3702 bp::list python_normals;
3704 for (
int idx = 0; idx < 3; ++idx) {
3705 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3706 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3709 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3710 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3712 std::regex reg_disp_name(disp_block_name);
3713 std::regex reg_traction_name(traction_block_name);
3715 np::ndarray np_analytical_expr = np::empty(
3716 bp::make_tuple(nb_gauss_pts, 3), np::dtype::get_builtin<double>());
3718 if (std::regex_match(block_name, reg_disp_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_disp() python call");
3728 }
else if (std::regex_match(block_name, reg_traction_name)) {
3730 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3731 bp::extract<np::ndarray>(python_coords[1]),
3732 bp::extract<np::ndarray>(python_coords[2]),
3733 bp::extract<np::ndarray>(python_normals[0]),
3734 bp::extract<np::ndarray>(python_normals[1]),
3735 bp::extract<np::ndarray>(python_normals[2]),
3736 block_name, np_analytical_expr),
3737 "Failed analytical_traction() python call");
3740 "Unknown analytical block");
3744 if (PetscUnlikely(np_analytical_expr.get_shape()[0] != nb_gauss_pts ||
3745 np_analytical_expr.get_shape()[1] != 3)) {
3748 "Wrong shape of analytical expression returned from "
3749 "python, expected: (" +
3750 std::to_string(nb_gauss_pts) +
", 3), got: (" +
3751 std::to_string(np_analytical_expr.get_shape()[0]) +
", " +
3752 std::to_string(np_analytical_expr.get_shape()[1]) +
")");
3755 double *analytical_expr_val_ptr =
3756 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3759 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3760 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3761 for (
int idx = 0; idx < 3; ++idx)
3762 v_analytical_expr(idx, gg) =
3763 *(analytical_expr_val_ptr + (3 * gg + idx));
3765 return v_analytical_expr;
3776 const std::string block_name) {
3778#ifdef ENABLE_PYTHON_BINDING
3779 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3783 bp::list python_coords;
3785 for (
int idx = 0; idx < 3; ++idx) {
3786 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3789 auto externalstrain_block_name =
"(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
3791 std::regex reg_externalstrain_name(externalstrain_block_name);
3793 np::ndarray np_analytical_expr = np::empty(
3794 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
3796 if (std::regex_match(block_name, reg_externalstrain_name)) {
3797 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalExternalStrain(
3798 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3799 bp::extract<np::ndarray>(python_coords[1]),
3800 bp::extract<np::ndarray>(python_coords[2]), block_name,
3801 np_analytical_expr),
3802 "Failed analytical_external_strain() python call");
3805 "Unknown analytical block");
3809 if (PetscUnlikely(np_analytical_expr.get_shape()[0] != nb_gauss_pts ||
3810 np_analytical_expr.get_shape()[1] != 1)) {
3813 "Wrong shape of analytical expression returned from "
3814 "python, expected: (" +
3815 std::to_string(nb_gauss_pts) +
", 1), got: (" +
3816 std::to_string(np_analytical_expr.get_shape()[0]) +
", " +
3817 std::to_string(np_analytical_expr.get_shape()[1]) +
")");
3820 double *analytical_expr_val_ptr =
3821 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3824 v_analytical_expr.resize(nb_gauss_pts,
false);
3825 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg)
3826 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
3828 return v_analytical_expr;
3837 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3838 boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
3839 boost::shared_ptr<Range> ents_ptr)
3840 :
OP(broken_base_side_data, ents_ptr) {
3841 this->sourceVec = vec;
3842 this->betaCoeff = beta_coeff;
3846 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3847 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3848 :
OP(broken_base_side_data, ents_ptr) {
3849 this->sourceVec = boost::shared_ptr<MatrixDouble>();
3850 this->betaCoeff = beta_coeff;
3859 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3864 if (!brokenBaseSideData) {
3869 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
3877 OP::nbIntegrationPts = OP::getGaussPts().size2();
3879 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
3881 OP::locF.resize(OP::nbRows,
false);
3884 CHKERR this->iNtegrate(row_data);
3886 CHKERR this->aSsemble(row_data);
3890 switch (OP::opType) {
3892 for (
auto &bd : *brokenBaseSideData) {
3894 boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
3895 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
3896 this->sourceVec.reset();
3901 (std::string(
"wrong op type ") +
3902 OpBaseDerivativesBase::OpTypeNames[OP::opType])
3910 const std::string row_field,
3911 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3912 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3914 brokenBaseSideDataPtr(broken_base_side_data) {
3915 this->betaCoeff = beta_coeff;
3921 for (
auto &bd : (*brokenBaseSideDataPtr)) {
3926 if (this->sourceVec->size1() !=
SPACE_DIM) {
3928 "Inconsistent size of the source vector");
3930 if (this->sourceVec->size2() != OP::getGaussPts().size2()) {
3932 "Inconsistent size of the source vector");
3936 CHKERR OP::iNtegrate(data);
3938 this->sourceVec.reset();
3944 std::string row_field,
3945 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3946 ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
3947 boost::shared_ptr<Range> ents_ptr)
3948 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
3950 this->betaCoeff = beta;
3955 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3956 ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
3957 :
OP(broken_base_side_data, ents_ptr) {
3959 this->betaCoeff = beta;
3960 OP::assembleTranspose =
false;
3961 OP::onlyTranspose =
false;
3970 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3975 if (!brokenBaseSideData) {
3980 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
3981 EntityType col_type,
3986 auto check_if_assemble_transpose = [&] {
3988 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
3992 }
else if (OP::assembleTranspose) {
3998 OP::rowSide = row_side;
3999 OP::rowType = row_type;
4000 OP::colSide = col_side;
4001 OP::colType = col_type;
4003 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
4005 CHKERR this->iNtegrate(row_data, col_data);
4006 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
4010 switch (OP::opType) {
4013 for (
auto &bd : *brokenBaseSideData) {
4016 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
4018 "base functions not set");
4022 OP::nbRows = bd.getData().getIndices().size();
4025 OP::nbIntegrationPts = OP::getGaussPts().size2();
4026 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
4034 bd.getSide(), bd.getSide(),
4037 bd.getType(), bd.getType(),
4040 bd.getData(), bd.getData()
4049 (std::string(
"wrong op type ") +
4050 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)
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)
