13#include <boost/math/constants/constants.hpp>
22#ifdef ENABLE_PYTHON_BINDING
23#include <boost/python.hpp>
24#include <boost/python/def.hpp>
25#include <boost/python/numpy.hpp>
26namespace bp = boost::python;
27namespace np = boost::python::numpy;
39 int nb_integration_pts = getGaussPts().size2();
41 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
42 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
47 auto t_eshelby_stress =
48 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
52 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
53 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
69 int nb_integration_pts = getGaussPts().size2();
85 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
86 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
87 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
89 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
90 dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts,
false);
91 dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 *
size_symm,
92 nb_integration_pts,
false);
93 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
95 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
96 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
97 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
98 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
100 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
101 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
102 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
103 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
104 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
105 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
106 dataAtPts->eigenValsC.resize(3, nb_integration_pts,
false);
107 dataAtPts->eigenVecsC.resize(9, nb_integration_pts,
false);
108 dataAtPts->nbUniqC.resize(nb_integration_pts,
false);
110 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
111 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
113 dataAtPts->internalStressAtPts.resize(9, nb_integration_pts,
false);
114 dataAtPts->internalStressAtPts.clear();
117 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
118 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
120 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
121 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
122 auto t_levi_kirchhoff_domega =
123 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
124 auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
125 dataAtPts->leviKirchhoffdLogStreatchAtPts);
126 auto t_levi_kirchhoff_dP =
127 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
128 auto t_approx_P_adjoint__dstretch =
129 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
130 auto t_approx_P_adjoint_log_du =
131 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
132 auto t_approx_P_adjoint_log_du_dP =
133 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
134 auto t_approx_P_adjoint_log_du_domega =
135 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
136 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
137 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
139 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
140 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
141 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
142 auto &nbUniq = dataAtPts->nbUniq;
145 auto t_eigen_vals_C = getFTensor1FromMat<3>(dataAtPts->eigenValsC);
146 auto t_eigen_vecs_C = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecsC);
147 auto &nbUniqC = dataAtPts->nbUniqC;
151 auto t_log_stretch_total =
152 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
154 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
157 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
158 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
159 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
161 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
169 ++t_levi_kirchhoff_domega;
170 ++t_levi_kirchhoff_dstreach;
171 ++t_levi_kirchhoff_dP;
172 ++t_approx_P_adjoint__dstretch;
173 ++t_approx_P_adjoint_log_du;
174 ++t_approx_P_adjoint_log_du_dP;
175 ++t_approx_P_adjoint_log_du_domega;
186 ++t_log_stretch_total;
197 auto bound_eig = [&](
auto &eig) {
199 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
200 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
201 for (
int ii = 0; ii != 3; ++ii) {
202 eig(ii) = std::min(std::max(zero, eig(ii)), large);
207 auto calculate_log_stretch = [&]() {
210 eigen_vec(
i,
j) = t_log_u(
i,
j);
212 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
216 t_nb_uniq = get_uniq_nb<3>(&eig(0));
220 t_eigen_vals(
i) = eig(
i);
221 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
224 auto get_t_diff_u = [&]() {
229 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
231 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
235 auto calculate_total_stretch = [&](
auto &t_h1) {
238 t_log_u2_h1(
i,
j) = 0;
239 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
244 t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
246 return std::make_pair(t_R_h1, t_inv_u_h1);
254 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
255 eigen_vec(
i,
j) = t_C_h1(
i,
j);
257 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
260 t_nb_uniq_C = get_uniq_nb<3>(&eig(0));
261 if (t_nb_uniq_C < 3) {
267 t_eigen_vals_C(
i) = eig(
i);
268 t_eigen_vecs_C(
i,
j) = eigen_vec(
i,
j);
272 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
274 auto f_inv_sqrt = [](
auto e) {
return 1. / std::sqrt(e); };
278 t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
280 return std::make_pair(t_R_h1, t_inv_u_h1);
284 auto no_h1_loop = [&]() {
294 "rotSelector should be large");
297 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
305 auto t_Ldiff_u = calculate_log_stretch();
308 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
324 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
325 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
327 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
328 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
333 "rotationSelector not handled");
336 constexpr double half_r = 1 / 2.;
337 constexpr double half_l = 1 / 2.;
340 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
343 t_approx_P_adjoint__dstretch(
l,
k) =
344 (t_R(
i,
l) * t_approx_P(
i,
k) + t_R(
i,
k) * t_approx_P(
i,
l));
345 t_approx_P_adjoint__dstretch(
l,
k) /= 2.;
347 t_approx_P_adjoint_log_du(
L) =
348 (t_approx_P_adjoint__dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L) +
349 t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
l,
k,
L) +
350 t_approx_P_adjoint__dstretch(
l,
k) * t_Ldiff_u(
k,
l,
L) +
351 t_approx_P_adjoint__dstretch(
k,
l) * t_Ldiff_u(
k,
l,
L)) /
355 t_levi_kirchhoff(
m) =
357 half_r * (t_diff_Rr(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
358 t_diff_Rr(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)))
362 half_l * (t_diff_Rl(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
363 t_diff_Rl(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)));
364 t_levi_kirchhoff(
m) /= 2.;
369 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
373 half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
375 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
378 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
380 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
381 t_R(
i,
l) * (t_Ldiff_u(
l,
k,
L) + t_Ldiff_u(
k,
l,
L)) / 2.;
385 t_A(
k,
l,
m) = t_diff_Rr(
i,
l,
m) * t_approx_P(
i,
k) +
386 t_diff_Rr(
i,
k,
m) * t_approx_P(
i,
l);
389 t_B(
k,
l,
m) = t_diff_Rl(
i,
l,
m) * t_approx_P(
i,
k) +
390 t_diff_Rl(
i,
k,
m) * t_approx_P(
i,
l);
393 t_approx_P_adjoint_log_du_domega(
m,
L) =
394 half_r * (t_A(
k,
l,
m) *
395 (t_Ldiff_u(
k,
l,
L) + t_Ldiff_u(
k,
l,
L)) / 2) +
396 half_l * (t_B(
k,
l,
m) *
397 (t_Ldiff_u(
k,
l,
L) + t_Ldiff_u(
k,
l,
L)) / 2);
400 t_A(
k,
l,
m) = t_diff_R(
i,
l,
m) * t_approx_P(
i,
k);
401 t_approx_P_adjoint_log_du_domega(
m,
L) =
402 t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L);
405 t_levi_kirchhoff_dstreach(
m,
L) =
407 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)))
412 (t_diff_Rl(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k)));
414 t_levi_kirchhoff_dP(
m,
i,
k) =
415 half_r * t_diff_Rr(
i,
l,
m) * (t_u(
l,
k))
419 half_l * t_diff_Rl(
i,
l,
m) * (t_u(
l,
k));
421 t_levi_kirchhoff_domega(
m,
n) =
423 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
424 t_diff_diff_Rr(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)))
429 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
430 t_diff_diff_Rl(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)));
431 t_levi_kirchhoff_domega(
m,
n) /= 2.;
440 auto large_loop = [&]() {
450 "rotSelector should be large");
453 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
467 "gradApproximator not handled");
471 auto t_Ldiff_u = calculate_log_stretch();
473 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
476 t_h_u(
l,
k) = t_u(
l, o) * t_h1(o,
k);
478 t_Ldiff_h_u(
l,
k,
L) = t_Ldiff_u(
l, o,
L) * t_h1(o,
k);
494 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
495 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
497 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
498 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
501 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
502 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
504 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
505 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
507 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
508 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
513 "rotationSelector not handled");
516 constexpr double half_r = 1 / 2.;
517 constexpr double half_l = 1 / 2.;
520 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
523 t_approx_P_adjoint__dstretch(
l, o) =
524 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
525 t_approx_P_adjoint_log_du(
L) =
526 t_approx_P_adjoint__dstretch(
l, o) * t_Ldiff_u(
l, o,
L);
529 t_levi_kirchhoff(
m) =
531 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
535 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
540 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
544 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
546 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
549 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
551 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
552 t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
556 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
558 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
560 t_approx_P_adjoint_log_du_domega(
m,
L) =
561 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
565 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
568 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
569 t_approx_P_adjoint_log_du_domega(
m,
L) =
570 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
573 t_levi_kirchhoff_dstreach(
m,
L) =
575 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
579 half_l * (t_diff_Rl(
i,
l,
m) *
580 (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
582 t_levi_kirchhoff_dP(
m,
i,
k) =
584 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
588 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
590 t_levi_kirchhoff_domega(
m,
n) =
592 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
597 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
606 auto moderate_loop = [&]() {
616 "rotSelector should be large");
619 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
630 "gradApproximator not handled");
634 auto t_Ldiff_u = calculate_log_stretch();
636 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
639 t_h_u(
l,
k) = (t_symm_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
641 t_L_h(
l,
k,
L) = t_L(
l, o,
L) * t_h1(o,
k);
657 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
658 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
660 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
661 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
664 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
665 t_diff_R(
i,
j,
l) = levi_civita(
i,
j,
l);
667 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
668 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
670 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
671 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
676 "rotationSelector not handled");
679 constexpr double half_r = 1 / 2.;
680 constexpr double half_l = 1 / 2.;
683 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
686 t_approx_P_adjoint__dstretch(
l, o) =
687 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
689 t_approx_P_adjoint_log_du(
L) =
690 t_approx_P_adjoint__dstretch(
l, o) * t_L(
l, o,
L);
693 t_levi_kirchhoff(
m) =
695 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
699 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
704 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
708 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
710 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
713 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
715 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
719 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_L_h(
l,
k,
L);
721 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_L_h(
l,
k,
L);
722 t_approx_P_adjoint_log_du_domega(
m,
L) =
723 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
727 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
730 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_L_h(
l,
k,
L);
731 t_approx_P_adjoint_log_du_domega(
m,
L) =
732 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
735 t_levi_kirchhoff_dstreach(
m,
L) =
736 half_r * (t_diff_Rr(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)))
740 half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
742 t_levi_kirchhoff_dP(
m,
i,
k) =
744 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
748 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
750 t_levi_kirchhoff_domega(
m,
n) =
752 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
757 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
766 auto small_loop = [&]() {
773 "rotSelector should be small");
776 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
785 "gradApproximator not handled");
792 t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
794 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
795 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
797 t_log_u2_h1(
i,
j) = 0;
798 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
801 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
803 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
804 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
807 t_approx_P_adjoint__dstretch(
i,
j) = t_approx_P(
i,
j);
808 t_approx_P_adjoint_log_du(
L) =
809 t_approx_P_adjoint__dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
810 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
811 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
814 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
815 t_levi_kirchhoff_dstreach(
m,
L) = 0;
816 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
817 t_levi_kirchhoff_domega(
m,
n) = 0;
843 "gradApproximator not handled");
856 auto N_InLoop = getNinTheLoop();
857 auto sensee = getSkeletonSense();
858 auto nb_gauss_pts = getGaussPts().size2();
859 auto t_normal = getFTensor1NormalsAtGaussPts();
862 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
864 dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
865 dataAtPts->tractionAtPts.clear();
867 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
869 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
870 t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
882 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
888 int nb_integration_pts = getGaussPts().size2();
889 auto t_w = getFTensor0IntegrationWeight();
890 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
891 auto t_coords = getFTensor1CoordsAtGaussPts();
892 auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
900 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
901 loc_reaction_forces(
i) += (t_traction(
i)) * t_w * getMeasure();
902 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
904 loc_moment_forces(
i) +=
905 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
906 t_traction(
k) * t_w * getMeasure();
913 reactionVec[0] += loc_reaction_forces(0);
914 reactionVec[1] += loc_reaction_forces(1);
915 reactionVec[2] += loc_reaction_forces(2);
916 reactionVec[3] += loc_moment_forces(0);
917 reactionVec[4] += loc_moment_forces(1);
918 reactionVec[5] += loc_moment_forces(2);
926 int nb_integration_pts = data.
getN().size1();
927 auto v = getVolume();
928 auto t_w = getFTensor0IntegrationWeight();
929 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
930 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
931 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
932 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
933 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
934 dataAtPts->wL2DotDotAtPts.clear();
936 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
938 auto piola_scale = dataAtPts->piolaScale;
939 auto alpha_w = alphaW / piola_scale;
940 auto alpha_rho = alphaRho / piola_scale;
942 int nb_base_functions = data.
getN().size2();
946 auto get_ftensor1 = [](
auto &
v) {
951 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
953 auto t_nf = get_ftensor1(nF);
955 for (; bb != nb_dofs / 3; ++bb) {
956 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
957 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
958 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
962 for (; bb != nb_base_functions; ++bb)
976 int nb_integration_pts = getGaussPts().size2();
977 auto v = getVolume();
978 auto t_w = getFTensor0IntegrationWeight();
979 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
980 auto t_omega_grad_dot =
981 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
982 int nb_base_functions = data.
getN().size2();
988 auto get_ftensor1 = [](
auto &
v) {
994 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
997 auto t_nf = get_ftensor1(nF);
999 for (; bb != nb_dofs / 3; ++bb) {
1000 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
1001 t_nf(
k) += (
a * alphaOmega ) *
1002 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
1007 for (; bb != nb_base_functions; ++bb) {
1021 int nb_integration_pts = data.
getN().size1();
1022 auto v = getVolume();
1023 auto t_w = getFTensor0IntegrationWeight();
1025 int nb_base_functions = data.
getN().size2() / 3;
1033 auto get_ftensor1 = [](
auto &
v) {
1040 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1041 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1043 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1045 auto t_nf = get_ftensor1(nF);
1050 for (; bb != nb_dofs / 3; ++bb) {
1051 t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
1052 t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
1053 t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
1058 for (; bb != nb_base_functions; ++bb)
1068 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1070 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1072 auto t_nf = get_ftensor1(nF);
1080 for (; bb != nb_dofs / 3; ++bb) {
1081 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1086 for (; bb != nb_base_functions; ++bb)
1100 int nb_integration_pts = data.
getN().size1();
1101 auto v = getVolume();
1102 auto t_w = getFTensor0IntegrationWeight();
1104 int nb_base_functions = data.
getN().size2() / 9;
1112 auto get_ftensor0 = [](
auto &
v) {
1118 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1119 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1121 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1123 auto t_nf = get_ftensor0(nF);
1128 for (; bb != nb_dofs; ++bb) {
1129 t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
1130 t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
1134 for (; bb != nb_base_functions; ++bb) {
1144 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1146 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1148 auto t_nf = get_ftensor0(nF);
1152 t_residuum(
i,
j) = t_h(
i,
j);
1155 for (; bb != nb_dofs; ++bb) {
1156 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1160 for (; bb != nb_base_functions; ++bb) {
1174 int nb_integration_pts = data.
getN().size1();
1175 auto v = getVolume();
1176 auto t_w = getFTensor0IntegrationWeight();
1177 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1178 int nb_base_functions = data.
getN().size2() / 3;
1181 auto get_ftensor1 = [](
auto &
v) {
1186 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1188 auto t_nf = get_ftensor1(nF);
1190 for (; bb != nb_dofs / 3; ++bb) {
1191 double div_row_base = t_row_diff_base_fun(
i,
i);
1192 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1194 ++t_row_diff_base_fun;
1196 for (; bb != nb_base_functions; ++bb) {
1197 ++t_row_diff_base_fun;
1211 int nb_integration_pts = getGaussPts().size2();
1214 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1216 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1217 if (tag_length != 9) {
1219 "Number of internal stress components should be 9 but is %d",
1224 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1225 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1226 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1227 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1229 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1230 dataAtPts->internalStressAtPts.clear();
1231 auto t_internal_stress =
1232 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1235 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1236 t_internal_stress(
L) = t_const_stress(
L);
1237 ++t_internal_stress;
1249 int nb_integration_pts = getGaussPts().size2();
1252 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1254 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1255 if (tag_length != 9) {
1257 "Number of internal stress components should be 9 but is %d",
1261 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1264 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1267 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1270 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1271 dataAtPts->internalStressAtPts.clear();
1272 auto t_internal_stress =
1273 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1278 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1279 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1280 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1281 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1285 ++t_internal_stress;
1296 int nb_integration_pts = data.
getN().size1();
1297 auto v = getVolume();
1298 auto t_w = getFTensor0IntegrationWeight();
1303 auto get_ftensor2 = [](
auto &
v) {
1305 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1308 auto t_internal_stress =
1309 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1314 int nb_base_functions = data.
getN().size2();
1316 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1318 auto t_nf = get_ftensor2(nF);
1321 t_symm_stress(
i,
j) =
1322 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1325 t_residual(
L) = t_L(
i,
j,
L) * t_symm_stress(
i,
j);
1328 for (; bb != nb_dofs / 6; ++bb) {
1329 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1333 for (; bb != nb_base_functions; ++bb)
1337 ++t_internal_stress;
1347 int nb_integration_pts = data.
getN().size1();
1348 auto v = getVolume();
1349 auto t_w = getFTensor0IntegrationWeight();
1351 auto get_ftensor2 = [](
auto &
v) {
1353 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1356 auto t_internal_stress =
1357 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1364 int nb_base_functions = data.
getN().size2();
1366 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1368 auto t_nf = get_ftensor2(nF);
1371 t_residual(
L) = t_L(M,
L) * t_internal_stress(M);
1374 for (; bb != nb_dofs / 6; ++bb) {
1375 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1379 for (; bb != nb_base_functions; ++bb)
1383 ++t_internal_stress;
1388template <AssemblyType A>
1394 for (
auto &bc : (*bcDispPtr)) {
1396 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1399 int nb_integration_pts = OP::getGaussPts().size2();
1400 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1401 auto t_w = OP::getFTensor0IntegrationWeight();
1402 int nb_base_functions = data.
getN().size2() / 3;
1409 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1411 scale *= scalingMethodsMap.at(bc.blockName)
1414 scale *= scalingMethodsMap.at(bc.blockName)
1415 ->getScale(OP::getFEMethod()->ts_t);
1419 <<
"No scaling method found for " << bc.blockName;
1426 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1427 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1429 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1431 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1435 for (; bb != nb_base_functions; ++bb)
1447 return OP::iNtegrate(data);
1450template <AssemblyType A>
1458 double time = OP::getFEMethod()->ts_t;
1466 for (
auto &bc : (*bcRotPtr)) {
1468 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1470 int nb_integration_pts = OP::getGaussPts().size2();
1471 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1472 auto t_w = OP::getFTensor0IntegrationWeight();
1474 int nb_base_functions = data.
getN().size2() / 3;
1477 auto get_ftensor1 = [](
auto &
v) {
1490 auto get_rotation_angle = [&]() {
1491 double theta = bc.theta;
1492 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1493 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1498 auto get_rotation = [&](
auto theta) {
1500 if (bc.vals.size() == 7) {
1501 t_omega(0) = bc.vals[4];
1502 t_omega(1) = bc.vals[5];
1503 t_omega(2) = bc.vals[6];
1506 t_omega(
i) = OP::getFTensor1Normal()(
i);
1509 t_omega(
i) *= theta;
1516 auto t_R = get_rotation(get_rotation_angle());
1517 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1519 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1521 t_delta(
i) = t_center(
i) - t_coords(
i);
1523 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1525 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1527 for (; bb != nb_dofs / 3; ++bb) {
1528 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1532 for (; bb != nb_base_functions; ++bb)
1545 return OP::iNtegrate(data);
1548template <AssemblyType A>
1552 double time = OP::getFEMethod()->ts_t;
1560 for (
auto &bc : (*bcDispPtr)) {
1562 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1564 for (
auto &bd : (*brokenBaseSideDataPtr)) {
1566 auto t_approx_P = getFTensor2FromMat<3, 3>(bd.getFlux());
1567 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1568 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1569 auto t_w = OP::getFTensor0IntegrationWeight();
1577 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1578 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1581 <<
"No scaling method found for " << bc.blockName;
1585 double val =
scale * bc.val;
1588 int nb_integration_pts = OP::getGaussPts().size2();
1589 int nb_base_functions = data.
getN().size2();
1591 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1594 t_N(
i) = t_normal(
i);
1598 t_P(
i,
j) = t_N(
i) * t_N(
j);
1603 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
1607 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
1609 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1611 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1612 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
1616 for (; bb != nb_base_functions; ++bb)
1630template <AssemblyType A>
1636 double time = OP::getFEMethod()->ts_t;
1641 int row_nb_dofs = row_data.
getIndices().size();
1642 int col_nb_dofs = col_data.
getIndices().size();
1643 auto &locMat = OP::locMat;
1644 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1650 for (
auto &bc : (*bcDispPtr)) {
1652 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1654 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1655 auto t_w = OP::getFTensor0IntegrationWeight();
1661 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1662 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1665 <<
"No scaling method found for " << bc.blockName;
1668 int nb_integration_pts = OP::getGaussPts().size2();
1669 int row_nb_dofs = row_data.
getIndices().size();
1670 int col_nb_dofs = col_data.
getIndices().size();
1671 int nb_base_functions = row_data.
getN().size2();
1676 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1679 t_N(
i) = t_normal(
i);
1683 t_P(
i,
j) = t_N(
i) * t_N(
j);
1686 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1689 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1690 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1693 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1694 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1701 for (; rr != nb_base_functions; ++rr)
1708 locMat *= OP::getMeasure();
1715template <AssemblyType A>
1721 double time = OP::getFEMethod()->ts_t;
1726 int row_nb_dofs = row_data.
getIndices().size();
1727 int col_nb_dofs = col_data.
getIndices().size();
1728 auto &locMat = OP::locMat;
1729 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1735 for (
auto &bc : (*bcDispPtr)) {
1737 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1739 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1740 auto t_w = OP::getFTensor0IntegrationWeight();
1749 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1750 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1753 <<
"No scaling method found for " << bc.blockName;
1756 int nb_integration_pts = OP::getGaussPts().size2();
1757 int nb_base_functions = row_data.
getN().size2();
1760 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1763 t_N(
i) = t_normal(
i);
1767 t_P(
i,
j) = t_N(
i) * t_N(
j);
1772 t_d_res(
i,
j) = t_Q(
i,
j);
1775 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1776 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1779 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1781 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1788 for (; rr != nb_base_functions; ++rr)
1795 locMat *= OP::getMeasure();
1802 return OP::iNtegrate(data);
1807 return OP::iNtegrate(row_data, col_data);
1812 return OP::iNtegrate(row_data, col_data);
1815template <AssemblyType A>
1819 double time = OP::getFEMethod()->ts_t;
1827 for (
auto &bc : (*bcDispPtr)) {
1829 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1834 auto [block_name, v_analytical_expr] =
1839 int nb_integration_pts = OP::getGaussPts().size2();
1840 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1841 auto t_w = OP::getFTensor0IntegrationWeight();
1842 int nb_base_functions = data.
getN().size2() / 3;
1844 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1850 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1852 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1853 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1856 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1858 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1862 for (; bb != nb_base_functions; ++bb)
1876 return OP::iNtegrate(data);
1885 int nb_integration_pts = getGaussPts().size2();
1886 int nb_base_functions = data.
getN().size2();
1888 double time = getFEMethod()->ts_t;
1894 if (this->locF.size() != nb_dofs)
1896 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1899 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1902 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1904 auto t_w = getFTensor0IntegrationWeight();
1905 auto t_coords = getFTensor1CoordsAtGaussPts();
1907 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1909 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1911 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1912 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1914 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1915 t_f(
i) -= (time_scale * t_w * t_row_base * tau) * (t_val(
i) *
scale);
1920 for (; rr != nb_base_functions; ++rr)
1925 this->locF *= getMeasure();
1931 for (
auto &bc : *(bcData)) {
1932 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1934 double time_scale = 1;
1935 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1936 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1942 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
1946 return -y * (y - 1) / 0.25;
1948 CHKERR integrate_rhs(bc, calc_tau, time_scale);
1951 bc, [](
double,
double,
double) {
return 1; }, time_scale);
1965 int nb_integration_pts = getGaussPts().size2();
1966 int nb_base_functions = data.
getN().size2();
1968 double time = getFEMethod()->ts_t;
1974 if (this->locF.size() != nb_dofs)
1976 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1979 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1984 auto t_w = getFTensor0IntegrationWeight();
1985 auto t_coords = getFTensor1CoordsAtGaussPts();
1986 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
1987 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
1989 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
1991 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1993 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2005 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2008 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2009 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2010 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2012 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2014 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2016 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2018 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2019 t_f(
i) += t_row_base * t_val(
i);
2024 for (; rr != nb_base_functions; ++rr)
2039 for (
auto &bc : *(bcData)) {
2040 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2042 double time_scale = 1;
2043 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2044 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2050 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2057template <AssemblyType A>
2068 double time = OP::getFEMethod()->ts_t;
2073 int nb_base_functions = row_data.
getN().size2();
2074 int row_nb_dofs = row_data.
getIndices().size();
2075 int col_nb_dofs = col_data.
getIndices().size();
2076 int nb_integration_pts = OP::getGaussPts().size2();
2077 auto &locMat = OP::locMat;
2078 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2081auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2086 auto t_w = OP::getFTensor0IntegrationWeight();
2087 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2088 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2089 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2091 auto t_grad_gamma_u = getFTensor2FromMat<3, 2>(*hybridGradDispPtr);
2094 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2104 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2105 auto t_val = time_scale * t_w * tau * val;
2108 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2109 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2112 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2114 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2115 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2116 t_kd(
j,
l) * t_diff_col_base(N0)
2120 (FTensor::levi_civita<double>(
i,
j,
k) *
2121 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2122 t_kd(
k,
l) * t_diff_col_base(N1);
2124 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2131 for (; rr != nb_base_functions; ++rr)
2147 for (
auto &bc : *(bcData)) {
2148 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2150 double time_scale = 1;
2151 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2152 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2158 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2169 return OP::iNtegrate(row_data, col_data);
2179 int nb_integration_pts = getGaussPts().size2();
2180 int nb_base_functions = data.
getN().size2();
2182 double time = getFEMethod()->ts_t;
2188 if (this->locF.size() != nb_dofs)
2190 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2195 for (
auto &bc : *(bcData)) {
2196 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2200 auto [block_name, v_analytical_expr] =
2203 getFTensor1FromMat<3>(v_analytical_expr);
2205 auto t_w = getFTensor0IntegrationWeight();
2206 auto t_coords = getFTensor1CoordsAtGaussPts();
2208 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2210 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2212 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2214 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2215 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2220 for (; rr != nb_base_functions; ++rr)
2226 this->locF *= getMeasure();
2235 int nb_integration_pts = row_data.
getN().size1();
2236 int row_nb_dofs = row_data.
getIndices().size();
2237 int col_nb_dofs = col_data.
getIndices().size();
2238 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2240 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2243 auto v = getVolume();
2244 auto t_w = getFTensor0IntegrationWeight();
2245 int row_nb_base_functions = row_data.
getN().size2();
2247 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2250 for (; rr != row_nb_dofs / 3; ++rr) {
2252 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2253 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2254 double div_col_base = t_col_diff_base_fun(
i,
i);
2255 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2257 ++t_col_diff_base_fun;
2261 for (; rr != row_nb_base_functions; ++rr)
2272 if (alphaW < std::numeric_limits<double>::epsilon() &&
2273 alphaRho < std::numeric_limits<double>::epsilon())
2276 const int nb_integration_pts = row_data.
getN().size1();
2277 const int row_nb_dofs = row_data.
getIndices().size();
2278 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2280 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2286 auto v = getVolume();
2287 auto t_w = getFTensor0IntegrationWeight();
2289 auto piola_scale = dataAtPts->piolaScale;
2290 auto alpha_w = alphaW / piola_scale;
2291 auto alpha_rho = alphaRho / piola_scale;
2293 int row_nb_base_functions = row_data.
getN().size2();
2296 double ts_scale = alpha_w * getTSa();
2297 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2298 ts_scale += alpha_rho * getTSaa();
2300 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2301 double a =
v * t_w * ts_scale;
2304 for (; rr != row_nb_dofs / 3; ++rr) {
2307 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2308 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2309 const double b =
a * t_row_base_fun * t_col_base_fun;
2318 for (; rr != row_nb_base_functions; ++rr)
2337 int nb_integration_pts = row_data.
getN().size1();
2338 int row_nb_dofs = row_data.
getIndices().size();
2339 int col_nb_dofs = col_data.
getIndices().size();
2340 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2343 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2345 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2347 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2349 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2351 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2353 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2356 auto v = getVolume();
2357 auto t_w = getFTensor0IntegrationWeight();
2359 int row_nb_base_functions = row_data.
getN().size2();
2362 auto t_approx_P_adjoint_log_du_dP =
2363 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2365 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2368 for (; rr != row_nb_dofs / 6; ++rr) {
2371 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2373 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2375 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2383 for (; rr != row_nb_base_functions; ++rr)
2386 ++t_approx_P_adjoint_log_du_dP;
2402 int nb_integration_pts = row_data.
getN().size1();
2403 int row_nb_dofs = row_data.
getIndices().size();
2404 int col_nb_dofs = col_data.
getIndices().size();
2405 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2407 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2411 auto v = getVolume();
2412 auto t_w = getFTensor0IntegrationWeight();
2415 int row_nb_base_functions = row_data.
getN().size2();
2417 auto t_approx_P_adjoint_log_du_dP =
2418 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2420 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2423 for (; rr != row_nb_dofs / 6; ++rr) {
2424 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2425 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2426 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2428 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2435 for (; rr != row_nb_base_functions; ++rr)
2438 ++t_approx_P_adjoint_log_du_dP;
2450 int row_nb_dofs = row_data.
getIndices().size();
2451 int col_nb_dofs = col_data.
getIndices().size();
2452 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2455 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2457 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2459 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2461 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2463 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2465 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2475 auto v = getVolume();
2476 auto t_w = getFTensor0IntegrationWeight();
2477 auto t_approx_P_adjoint_log_du_domega =
2478 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2480 int row_nb_base_functions = row_data.
getN().size2();
2483 int nb_integration_pts = row_data.
getN().size1();
2485 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2489 for (; rr != row_nb_dofs / 6; ++rr) {
2491 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2492 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2493 double v =
a * t_row_base_fun * t_col_base_fun;
2494 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2501 for (; rr != row_nb_base_functions; ++rr)
2505 ++t_approx_P_adjoint_log_du_domega;
2514 int nb_integration_pts = getGaussPts().size2();
2515 int row_nb_dofs = row_data.
getIndices().size();
2516 int col_nb_dofs = col_data.
getIndices().size();
2517 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2521 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2522 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2524 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2525 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2527 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2528 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2536 auto v = getVolume();
2537 auto t_w = getFTensor0IntegrationWeight();
2538 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2539 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2540 int row_nb_base_functions = row_data.
getN().size2();
2542 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2545 for (; rr != row_nb_dofs / 3; ++rr) {
2546 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2547 const double b =
a * t_row_base_fun;
2549 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2550 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2556 for (; rr != row_nb_base_functions; ++rr) {
2560 ++t_levi_kirchhoff_du;
2576 int nb_integration_pts = getGaussPts().size2();
2577 int row_nb_dofs = row_data.
getIndices().size();
2578 int col_nb_dofs = col_data.
getIndices().size();
2579 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2582 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2584 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2586 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2591 auto v = getVolume();
2592 auto t_w = getFTensor0IntegrationWeight();
2594 int row_nb_base_functions = row_data.
getN().size2();
2596 auto t_levi_kirchhoff_dP =
2597 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2599 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2602 for (; rr != row_nb_dofs / 3; ++rr) {
2603 double b =
a * t_row_base_fun;
2605 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2606 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2607 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2613 for (; rr != row_nb_base_functions; ++rr) {
2618 ++t_levi_kirchhoff_dP;
2626 int nb_integration_pts = getGaussPts().size2();
2627 int row_nb_dofs = row_data.
getIndices().size();
2628 int col_nb_dofs = col_data.
getIndices().size();
2630 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2632 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2639 auto v = getVolume();
2640 auto t_w = getFTensor0IntegrationWeight();
2641 auto t_levi_kirchoff_dP =
2642 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2644 int row_nb_base_functions = row_data.
getN().size2();
2647 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2650 for (; rr != row_nb_dofs / 3; ++rr) {
2651 double b =
a * t_row_base_fun;
2652 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2653 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2654 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2655 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2662 for (; rr != row_nb_base_functions; ++rr) {
2666 ++t_levi_kirchoff_dP;
2674 int nb_integration_pts = getGaussPts().size2();
2675 int row_nb_dofs = row_data.
getIndices().size();
2676 int col_nb_dofs = col_data.
getIndices().size();
2677 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2680 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2682 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2684 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2697 auto v = getVolume();
2698 auto ts_a = getTSa();
2699 auto t_w = getFTensor0IntegrationWeight();
2700 auto t_levi_kirchhoff_domega =
2701 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2702 int row_nb_base_functions = row_data.
getN().size2();
2707 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2709 double c = (
a * alphaOmega) * (ts_a );
2712 for (; rr != row_nb_dofs / 3; ++rr) {
2713 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2714 const double b =
a * t_row_base_fun;
2717 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2718 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2719 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2727 for (; rr != row_nb_base_functions; ++rr) {
2732 ++t_levi_kirchhoff_domega;
2740 if (dataAtPts->matInvD.size1() ==
size_symm &&
2741 dataAtPts->matInvD.size2() ==
size_symm) {
2754 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2757 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2759 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2761 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2766 int nb_integration_pts = getGaussPts().size2();
2767 int row_nb_dofs = row_data.
getIndices().size();
2768 int col_nb_dofs = col_data.
getIndices().size();
2770 auto v = getVolume();
2771 auto t_w = getFTensor0IntegrationWeight();
2772 int row_nb_base_functions = row_data.
getN().size2() / 3;
2779 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2780 &*dataAtPts->matInvD.data().begin());
2783 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2787 for (; rr != row_nb_dofs / 3; ++rr) {
2789 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2790 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2791 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2799 for (; rr != row_nb_base_functions; ++rr)
2812 if (dataAtPts->matInvD.size1() ==
size_symm &&
2813 dataAtPts->matInvD.size2() ==
size_symm) {
2827 int nb_integration_pts = getGaussPts().size2();
2828 int row_nb_dofs = row_data.
getIndices().size();
2829 int col_nb_dofs = col_data.
getIndices().size();
2831 auto v = getVolume();
2832 auto t_w = getFTensor0IntegrationWeight();
2833 int row_nb_base_functions = row_data.
getN().size2() / 9;
2840 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2841 &*dataAtPts->matInvD.data().begin());
2844 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2848 for (; rr != row_nb_dofs; ++rr) {
2850 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2852 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2859 for (; rr != row_nb_base_functions; ++rr)
2870 if (dataAtPts->matInvD.size1() ==
size_symm &&
2871 dataAtPts->matInvD.size2() ==
size_symm) {
2885 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2888 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
2893 int nb_integration_pts = getGaussPts().size2();
2894 int row_nb_dofs = row_data.
getIndices().size();
2895 int col_nb_dofs = col_data.
getIndices().size();
2897 auto v = getVolume();
2898 auto t_w = getFTensor0IntegrationWeight();
2899 int row_nb_base_functions = row_data.
getN().size2() / 9;
2908 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2909 &*dataAtPts->matInvD.data().begin());
2912 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2915 auto t_m = get_ftensor1(
K, 0, 0);
2918 for (; rr != row_nb_dofs; ++rr) {
2920 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2921 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2929 for (; rr != row_nb_base_functions; ++rr)
2948 int nb_integration_pts = row_data.
getN().size1();
2949 int row_nb_dofs = row_data.
getIndices().size();
2950 int col_nb_dofs = col_data.
getIndices().size();
2952 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2955 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2957 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2959 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2964 auto v = getVolume();
2965 auto t_w = getFTensor0IntegrationWeight();
2966 int row_nb_base_functions = row_data.
getN().size2() / 3;
2969 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2971 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2975 for (; rr != row_nb_dofs / 3; ++rr) {
2978 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2981 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2982 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2983 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
2991 for (; rr != row_nb_base_functions; ++rr)
3011 int nb_integration_pts = row_data.
getN().size1();
3012 int row_nb_dofs = row_data.
getIndices().size();
3013 int col_nb_dofs = col_data.
getIndices().size();
3015 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3017 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3020 auto v = getVolume();
3021 auto t_w = getFTensor0IntegrationWeight();
3022 int row_nb_base_functions = row_data.
getN().size2() / 9;
3025 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
3026 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3030 for (; rr != row_nb_dofs; ++rr) {
3033 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3036 auto t_m = get_ftensor2(
K, rr, 0);
3037 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3038 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3046 for (; rr != row_nb_base_functions; ++rr)
3059 if (tagSense != getSkeletonSense())
3062 auto get_tag = [&](
auto name) {
3063 auto &mob = getPtrFE()->mField.get_moab();
3069 auto get_tag_value = [&](
auto &&tag,
int dim) {
3070 auto &mob = getPtrFE()->mField.get_moab();
3071 auto face = getSidePtrFE()->getFEEntityHandle();
3072 std::vector<double> value(dim);
3073 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3077 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3078 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3080 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3081 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3086 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3087 Tag th_detF = create_tag(
"detF", 1);
3088 Tag th_traction = create_tag(
"traction", 3);
3089 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3091 Tag th_energy = create_tag(
"Energy", 1);
3093 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3094 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3095 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3097 auto t_normal = getFTensor1NormalsAtGaussPts();
3098 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3108 if (dataAtPts->energyAtPts.size() == 0) {
3110 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3111 dataAtPts->energyAtPts.clear();
3120 auto set_float_precision = [](
const double x) {
3121 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3128 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3130 v = set_float_precision(
v);
3131 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3138 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3141 for (
auto &
a :
v.data())
3142 a = set_float_precision(
a);
3143 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3144 &*
v.data().begin());
3152 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3154 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3156 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3158 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3160 t_m(
i,
j) = t_d(
i,
j);
3161 for (
auto &
v :
m.data())
3162 v = set_float_precision(
v);
3163 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3164 &*
m.data().begin());
3168 const auto nb_gauss_pts = getGaussPts().size2();
3169 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3172 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3174 CHKERR save_vec_tag(th_traction, t_traction, gg);
3176 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3177 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3178 CHKERR save_scal_tag(th_energy, t_energy, gg);
3182 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3183 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3184 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3193 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3194 std::vector<FieldSpace> spaces, std::string geom_field_name,
3195 boost::shared_ptr<Range> crack_front_edges_ptr) {
3198 constexpr bool scale_l2 =
false;
3202 "Scale L2 Ainsworth Legendre base is not implemented");
3211 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3212 std::vector<FieldSpace> spaces, std::string geom_field_name,
3213 boost::shared_ptr<Range> crack_front_edges_ptr) {
3216 constexpr bool scale_l2 =
false;
3220 "Scale L2 Ainsworth Legendre base is not implemented");
3229 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3230 std::vector<FieldSpace> spaces, std::string geom_field_name,
3231 boost::shared_ptr<Range> crack_front_edges_ptr,
3232 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3233 boost::shared_ptr<MatrixDouble> inv_jac) {
3237 auto jac = boost::make_shared<MatrixDouble>();
3238 auto det = boost::make_shared<VectorDouble>();
3240 geom_field_name, jac));
3246 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3248 if (scale_l2_ainsworth_legendre_base) {
3256 boost::shared_ptr<MatrixDouble> jac,
3257 boost::shared_ptr<Range> edges_ptr)
3266 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3269 return OP::doWork(side, type, data);
3274 boost::shared_ptr<Range> edgesPtr;
3277 auto jac = boost::make_shared<MatrixDouble>();
3278 auto det = boost::make_shared<VectorDouble>();
3280 geom_field_name, jac,
3282 : boost::make_shared<
Range>()));
3352 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3353 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3354 if (getNinTheLoop() == 0) {
3355 dataAtPts->faceMaterialForceAtPts.clear();
3356 dataAtPts->normalPressureAtPts.clear();
3358 auto loop_size = getLoopSize();
3359 if (loop_size == 1) {
3360 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3361 auto pstatus = numebered_fe_ptr->getPStatus();
3362 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3369 auto t_normal = getFTensor1NormalsAtGaussPts();
3370 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3371 dataAtPts->faceMaterialForceAtPts);
3374 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3378 getFTensor1FromMat<SPACE_DIM>(dataAtPts->hybridDispAtPts);
3379 auto t_grad_u_gamma =
3380 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3382 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3383 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
3405 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3406 t_N(
I) = t_normal(
I);
3409 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3411 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3414 t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3417 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3419 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3425 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3428 t_N(
I) = t_normal(
I);
3433 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3436 t_grad_u(
i,
J) = t_grad_u_gamma(
i,
J) + (2 * t_R(
i,
K) * t_N(
K) -
3437 (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3440 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3443 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3446 t_p += t_N(
I) * (t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J))) / loop_size;
3454 "Grffith energy release "
3455 "selector not implemented");
3459 auto side_fe_ptr = getSidePtrFE();
3460 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3461 auto pstatus = side_fe_mi_ptr->getPStatus();
3463 auto owner = side_fe_mi_ptr->getOwnerProc();
3465 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3466 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3467 << getNinTheLoop() <<
" loop size " << getLoopSize();
3479 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3480 auto pstatus = fe_mi_ptr->getPStatus();
3482 auto owner = fe_mi_ptr->getOwnerProc();
3484 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3485 <<
" " << getPtrFE()->mField.get_comm_rank();
3493 double face_pressure = 0.;
3494 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3495 dataAtPts->faceMaterialForceAtPts);
3498 auto t_w = getFTensor0IntegrationWeight();
3499 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3500 t_face_T(
I) += t_w * t_T(
I);
3501 face_pressure += t_w * t_p;
3506 t_face_T(
I) *= getMeasure();
3507 face_pressure *= getMeasure();
3509 auto get_tag = [&](
auto name,
auto dim) {
3510 auto &moab = getPtrFE()->mField.get_moab();
3512 double def_val[] = {0., 0., 0.};
3513 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3514 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3519 auto set_tag = [&](
auto &&tag,
auto ptr) {
3520 auto &moab = getPtrFE()->mField.get_moab();
3521 auto face = getPtrFE()->getFEEntityHandle();
3522 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3525 set_tag(get_tag(
"MaterialForce", 3), &t_face_T(0));
3526 set_tag(get_tag(
"FacePressure", 1), &face_pressure);
3532template <
typename OP_PTR>
3534 const std::string block_name) {
3536 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
3538 auto ts_time = op_ptr->getTStime();
3539 auto ts_time_step = op_ptr->getTStimeStep();
3545 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
3546 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
3548 auto v_analytical_expr =
3550 m_ref_coords, m_ref_normals, block_name);
3553 if (v_analytical_expr.size2() != nb_gauss_pts)
3555 "Wrong number of integration pts");
3558 return std::make_tuple(block_name, v_analytical_expr);
3562#ifdef ENABLE_PYTHON_BINDING
3564 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3569 auto main_module = bp::import(
"__main__");
3570 mainNamespace = main_module.attr(
"__dict__");
3571 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3573 analyticalDispFun = mainNamespace[
"analytical_disp"];
3574 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3575 analyticalExternalStrainFun = mainNamespace[
"analytical_external_strain"];
3577 }
catch (bp::error_already_set
const &) {
3587 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3588 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3589 const std::string &block_name, np::ndarray &analytical_expr
3596 analytical_expr = bp::extract<np::ndarray>(
3597 analyticalDispFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3599 }
catch (bp::error_already_set
const &) {
3609 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3610 np::ndarray nx, np::ndarray ny, np::ndarray nz,
3611 const std::string& block_name, np::ndarray &analytical_expr
3618 analytical_expr = bp::extract<np::ndarray>(
3619 analyticalTractionFun(delta_t,
t, x, y, z, nx, ny, nz, block_name));
3621 }
catch (bp::error_already_set
const &) {
3629 MoFEMErrorCode AnalyticalExprPython::evalAnalyticalExternalStrain(
3631 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3632 const std::string& block_name, np::ndarray &analytical_expr
3639 analytical_expr = bp::extract<np::ndarray>(
3640 analyticalExternalStrainFun(delta_t,
t, x, y, z, block_name));
3642 }
catch (bp::error_already_set
const &) {
3650boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3652inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3654 auto dtype = np::dtype::get_builtin<double>();
3655 auto size = bp::make_tuple(nb_gauss_pts);
3656 auto stride = bp::make_tuple(3 *
sizeof(
double));
3657 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3666 const std::string block_name) {
3668#ifdef ENABLE_PYTHON_BINDING
3669 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3674 bp::list python_coords;
3675 bp::list python_normals;
3677 for (
int idx = 0; idx < 3; ++idx) {
3678 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3679 python_normals.append(convert_to_numpy(v_ref_normals, nb_gauss_pts, idx));
3682 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3683 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3685 std::regex reg_disp_name(disp_block_name);
3686 std::regex reg_traction_name(traction_block_name);
3688 np::ndarray np_analytical_expr = np::empty(
3689 bp::make_tuple(nb_gauss_pts, 3), np::dtype::get_builtin<double>());
3691 if (std::regex_match(block_name, reg_disp_name)) {
3693 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3694 bp::extract<np::ndarray>(python_coords[1]),
3695 bp::extract<np::ndarray>(python_coords[2]),
3696 bp::extract<np::ndarray>(python_normals[0]),
3697 bp::extract<np::ndarray>(python_normals[1]),
3698 bp::extract<np::ndarray>(python_normals[2]),
3699 block_name, np_analytical_expr),
3700 "Failed analytical_disp() python call");
3701 }
else if (std::regex_match(block_name, reg_traction_name)) {
3703 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3704 bp::extract<np::ndarray>(python_coords[1]),
3705 bp::extract<np::ndarray>(python_coords[2]),
3706 bp::extract<np::ndarray>(python_normals[0]),
3707 bp::extract<np::ndarray>(python_normals[1]),
3708 bp::extract<np::ndarray>(python_normals[2]),
3709 block_name, np_analytical_expr),
3710 "Failed analytical_traction() python call");
3713 "Unknown analytical block");
3716 double *analytical_expr_val_ptr =
3717 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3720 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3721 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3722 for (
int idx = 0; idx < 3; ++idx)
3723 v_analytical_expr(idx, gg) =
3724 *(analytical_expr_val_ptr + (3 * gg + idx));
3726 return v_analytical_expr;
3737 const std::string block_name) {
3739#ifdef ENABLE_PYTHON_BINDING
3740 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3744 bp::list python_coords;
3746 for (
int idx = 0; idx < 3; ++idx) {
3747 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3750 auto externalstrain_block_name =
"(.*)ANALYTICAL_EXTERNALSTRAIN(.*)";
3752 std::regex reg_externalstrain_name(externalstrain_block_name);
3754 np::ndarray np_analytical_expr = np::empty(
3755 bp::make_tuple(nb_gauss_pts), np::dtype::get_builtin<double>());
3757 if (std::regex_match(block_name, reg_externalstrain_name)) {
3758 CHK_MOAB_THROW(analytical_expr_ptr->evalAnalyticalExternalStrain(
3759 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3760 bp::extract<np::ndarray>(python_coords[1]),
3761 bp::extract<np::ndarray>(python_coords[2]), block_name,
3762 np_analytical_expr),
3763 "Failed analytical_external_strain() python call");
3766 "Unknown analytical block");
3769 double *analytical_expr_val_ptr =
3770 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3773 v_analytical_expr.resize(nb_gauss_pts,
false);
3774 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg)
3775 v_analytical_expr[gg] = *(analytical_expr_val_ptr + gg);
3777 return v_analytical_expr;
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
constexpr int SPACE_DIM
[Define dimension]
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)
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 iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Caluclate face material force and normal pressure at gauss points.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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)
