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;
35 const std::string block_name);
37template <
typename OP_PTR>
39 const std::string block_name) {
41 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
43 auto ts_time = op_ptr->getTStime();
44 auto ts_time_step = op_ptr->getTStimeStep();
50 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
53 ts_time_step, ts_time, nb_gauss_pts, m_ref_coords, block_name);
56 if (v_analytical_expr.size2() != nb_gauss_pts)
58 "Wrong number of integration pts");
61 return std::make_tuple(block_name, v_analytical_expr);
71 int nb_integration_pts = getGaussPts().size2();
73 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
74 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
79 auto t_eshelby_stress =
80 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
84 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
85 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
101 int nb_integration_pts = getGaussPts().size2();
117 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
118 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
119 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
121 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
122 dataAtPts->leviKirchhoffdOmegaAtPts.resize(9, nb_integration_pts,
false);
123 dataAtPts->leviKirchhoffdLogStreatchAtPts.resize(3 *
size_symm,
124 nb_integration_pts,
false);
125 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
127 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
128 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
129 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
130 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
132 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
133 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
134 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
135 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
136 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
137 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
139 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
140 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
142 dataAtPts->internalStressAtPts.resize(9, nb_integration_pts,
false);
143 dataAtPts->internalStressAtPts.clear();
146 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
147 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
149 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
150 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
151 auto t_levi_kirchhoff_domega =
152 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
153 auto t_levi_kirchhoff_dstreach = getFTensor2FromMat<3, size_symm>(
154 dataAtPts->leviKirchhoffdLogStreatchAtPts);
155 auto t_levi_kirchhoff_dP =
156 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
157 auto t_approx_P_adjoint_dstretch =
158 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
159 auto t_approx_P_adjoint_log_du =
160 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
161 auto t_approx_P_adjoint_log_du_dP =
162 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
163 auto t_approx_P_adjoint_log_du_domega =
164 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
165 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
166 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
168 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
169 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
170 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
171 auto t_log_stretch_total =
172 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
174 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
175 auto &nbUniq = dataAtPts->nbUniq;
180 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
181 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
182 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
184 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
192 ++t_levi_kirchhoff_domega;
193 ++t_levi_kirchhoff_dstreach;
194 ++t_levi_kirchhoff_dP;
195 ++t_approx_P_adjoint_dstretch;
196 ++t_approx_P_adjoint_log_du;
197 ++t_approx_P_adjoint_log_du_dP;
198 ++t_approx_P_adjoint_log_du_domega;
206 ++t_log_stretch_total;
217 auto bound_eig = [&](
auto &eig) {
219 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
220 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
221 for (
int ii = 0; ii != 3; ++ii) {
222 eig(ii) = std::min(std::max(zero, eig(ii)), large);
227 auto calculate_log_stretch = [&]() {
230 eigen_vec(
i,
j) = t_log_u(
i,
j);
232 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
236 t_nb_uniq = get_uniq_nb<3>(&eig(0));
240 t_eigen_vals(
i) = eig(
i);
241 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
244 auto get_t_diff_u = [&]() {
249 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
251 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
255 auto calculate_total_stretch = [&](
auto &t_h1) {
258 t_log_u2_h1(
i,
j) = 0;
259 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
264 t_inv_u_h1(
i,
j) = t_symm_kd(
i,
j);
266 return std::make_pair(t_R_h1, t_inv_u_h1);
274 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
275 eigen_vec(
i,
j) = t_C_h1(
i,
j);
277 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
286 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
288 auto f_inv_sqrt = [](
auto e) {
return 1. / std::sqrt(e); };
292 t_R_h1(
i,
j) = t_h1(
i, o) * t_inv_u_h1(o,
j);
294 return std::make_pair(t_R_h1, t_inv_u_h1);
298 auto no_h1_loop = [&]() {
308 "rotSelector should be large");
311 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
319 auto t_Ldiff_u = calculate_log_stretch();
322 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
338 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
339 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
341 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
342 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
347 "rotationSelector not handled");
350 constexpr double half_r = 1 / 2.;
351 constexpr double half_l = 1 / 2.;
354 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
357 t_approx_P_adjoint_dstretch(
l,
k) =
358 (t_R(
i,
l) * t_approx_P(
i,
k) + t_R(
i,
k) * t_approx_P(
i,
l));
359 t_approx_P_adjoint_dstretch(
l,
k) /= 2.;
361 t_approx_P_adjoint_log_du(
L) =
362 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
365 t_levi_kirchhoff(
m) =
367 half_r * (t_diff_Rr(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
368 t_diff_Rr(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)))
372 half_l * (t_diff_Rl(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
373 t_diff_Rl(
i,
k,
m) * (t_u(
l,
k) * t_approx_P(
i,
l)));
374 t_levi_kirchhoff(
m) /= 2.;
379 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_u(
l,
k))
383 half_l * (t_diff_Rl(
i,
l,
m) * t_u(
l,
k));
386 "symmetrySelector not handled");
389 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
391 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
394 t_A(
k,
l,
m) = t_diff_Rr(
i,
l,
m) * t_approx_P(
i,
k) +
395 t_diff_Rr(
i,
k,
m) * t_approx_P(
i,
l);
398 t_B(
k,
l,
m) = t_diff_Rl(
i,
l,
m) * t_approx_P(
i,
k) +
399 t_diff_Rl(
i,
k,
m) * t_approx_P(
i,
l);
402 t_approx_P_adjoint_log_du_domega(
m,
L) =
403 half_r * (t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L)) +
404 half_l * (t_B(
k,
l,
m) * t_Ldiff_u(
k,
l,
L));
406 t_levi_kirchhoff_domega(
m,
n) =
408 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
409 t_diff_diff_Rr(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)))
414 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k)) +
415 t_diff_diff_Rl(
i,
k,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
l)));
416 t_levi_kirchhoff_domega(
m,
n) /= 2.;
425 auto large_loop = [&]() {
435 "rotSelector should be large");
438 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
452 "gradApproximator not handled");
456 auto t_Ldiff_u = calculate_log_stretch();
458 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
461 t_h_u(
l,
k) = t_u(
l, o) * t_h1(o,
k);
463 t_Ldiff_h_u(
l,
k,
L) = t_Ldiff_u(
l, o,
L) * t_h1(o,
k);
479 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
480 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
482 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
483 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
486 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
487 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
489 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
490 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
492 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
493 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
498 "rotationSelector not handled");
501 constexpr double half_r = 1 / 2.;
502 constexpr double half_l = 1 / 2.;
505 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
508 t_approx_P_adjoint_dstretch(
l, o) =
509 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
510 t_approx_P_adjoint_log_du(
L) =
511 t_approx_P_adjoint_dstretch(
l, o) * t_Ldiff_u(
l, o,
L);
514 t_levi_kirchhoff(
m) =
516 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
520 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
525 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
529 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
531 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
534 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
536 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
537 t_R(
i,
l) * t_Ldiff_h_u(
l,
k,
L);
541 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
543 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
545 t_approx_P_adjoint_log_du_domega(
m,
L) =
546 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
550 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
553 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_h_u(
l,
k,
L);
554 t_approx_P_adjoint_log_du_domega(
m,
L) =
555 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
558 t_levi_kirchhoff_dstreach(
m,
L) =
560 (t_diff_Rr(
i,
l,
m) * (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)))
564 half_l * (t_diff_Rl(
i,
l,
m) *
565 (t_Ldiff_h_u(
l,
k,
L) * t_approx_P(
i,
k)));
567 t_levi_kirchhoff_dP(
m,
i,
k) =
569 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
573 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
575 t_levi_kirchhoff_domega(
m,
n) =
577 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
582 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
591 auto moderate_loop = [&]() {
601 "rotSelector should be large");
604 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
615 "gradApproximator not handled");
619 auto t_Ldiff_u = calculate_log_stretch();
621 auto [t_R_h1, t_inv_u_h1] = calculate_total_stretch(t_h1);
624 t_h_u(
l,
k) = (t_symm_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
626 t_L_h(
l,
k,
L) = t_L(
l, o,
L) * t_h1(o,
k);
642 t_diff_Rr(
i,
j,
l) = t_R(
i,
k) * levi_civita(
k,
j,
l);
643 t_diff_diff_Rr(
i,
j,
l,
m) = t_diff_R(
i,
k,
m) * levi_civita(
k,
j,
l);
645 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
k,
l) * t_R(
k,
j);
646 t_diff_diff_Rl(
i,
j,
l,
m) = levi_civita(
i,
k,
l) * t_diff_R(
k,
j,
m);
649 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
650 t_diff_R(
i,
j,
l) = levi_civita(
i,
j,
l);
652 t_diff_Rr(
i,
j,
l) = levi_civita(
i,
j,
l);
653 t_diff_diff_Rr(
i,
j,
l,
m) = 0;
655 t_diff_Rl(
i,
j,
l) = levi_civita(
i,
j,
l);
656 t_diff_diff_Rl(
i,
j,
l,
m) = 0;
661 "rotationSelector not handled");
664 constexpr double half_r = 1 / 2.;
665 constexpr double half_l = 1 / 2.;
668 t_h(
i,
k) = t_R(
i,
l) * t_h_u(
l,
k);
671 t_approx_P_adjoint_dstretch(
l, o) =
672 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
674 t_approx_P_adjoint_log_du(
L) =
675 t_approx_P_adjoint_dstretch(
l, o) * t_L(
l, o,
L);
678 t_levi_kirchhoff(
m) =
680 half_r * ((t_diff_Rr(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)))
684 half_l * ((t_diff_Rl(
i,
l,
m) * (t_h_u(
l,
k)) * t_approx_P(
i,
k)));
689 t_h_domega(
i,
k,
m) = half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
693 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
695 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_h_u(
l,
k);
698 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
700 t_approx_P_adjoint_log_du_dP(
i,
k,
L) = t_R(
i,
l) * t_L_h(
l,
k,
L);
704 t_A(
m,
L,
i,
k) = t_diff_Rr(
i,
l,
m) * t_L_h(
l,
k,
L);
706 t_B(
m,
L,
i,
k) = t_diff_Rl(
i,
l,
m) * t_L_h(
l,
k,
L);
707 t_approx_P_adjoint_log_du_domega(
m,
L) =
708 half_r * (t_A(
m,
L,
i,
k) * t_approx_P(
i,
k))
712 half_l * (t_B(
m,
L,
i,
k) * t_approx_P(
i,
k));
715 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_L_h(
l,
k,
L);
716 t_approx_P_adjoint_log_du_domega(
m,
L) =
717 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
720 t_levi_kirchhoff_dstreach(
m,
L) =
721 half_r * (t_diff_Rr(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)))
725 half_l * (t_diff_Rl(
i,
l,
m) * (t_L_h(
l,
k,
L) * t_approx_P(
i,
k)));
727 t_levi_kirchhoff_dP(
m,
i,
k) =
729 half_r * (t_diff_Rr(
i,
l,
m) * t_h_u(
l,
k))
733 half_l * (t_diff_Rl(
i,
l,
m) * t_h_u(
l,
k));
735 t_levi_kirchhoff_domega(
m,
n) =
737 (t_diff_diff_Rr(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)))
742 (t_diff_diff_Rl(
i,
l,
m,
n) * (t_h_u(
l,
k) * t_approx_P(
i,
k)));
751 auto small_loop = [&]() {
758 "rotSelector should be small");
761 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
770 "gradApproximator not handled");
777 t_Ldiff_u(
i,
j,
L) = calculate_log_stretch()(
i,
j,
L);
779 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
780 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
782 t_log_u2_h1(
i,
j) = 0;
783 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
786 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
788 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
789 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
792 t_approx_P_adjoint_dstretch(
i,
j) = t_approx_P(
i,
j);
793 t_approx_P_adjoint_log_du(
L) =
794 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
795 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
796 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
799 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
800 t_levi_kirchhoff_dstreach(
m,
L) = 0;
801 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
802 t_levi_kirchhoff_domega(
m,
n) = 0;
828 "gradApproximator not handled");
841 auto N_InLoop = getNinTheLoop();
842 auto sensee = getSkeletonSense();
843 auto nb_gauss_pts = getGaussPts().size2();
844 auto t_normal = getFTensor1NormalsAtGaussPts();
847 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
849 dataAtPts->tractionAtPts.resize(
SPACE_DIM, nb_gauss_pts,
false);
850 dataAtPts->tractionAtPts.clear();
852 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
854 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
855 t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
867 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
873 int nb_integration_pts = getGaussPts().size2();
874 auto t_w = getFTensor0IntegrationWeight();
875 auto t_traction = getFTensor1FromMat<SPACE_DIM>(dataAtPts->tractionAtPts);
876 auto t_coords = getFTensor1CoordsAtGaussPts();
877 auto t_spatial_disp = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
885 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
886 loc_reaction_forces(
i) += (t_traction(
i)) * t_w * getMeasure();
887 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
889 loc_moment_forces(
i) +=
890 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
891 t_traction(
k) * t_w * getMeasure();
898 reactionVec[0] += loc_reaction_forces(0);
899 reactionVec[1] += loc_reaction_forces(1);
900 reactionVec[2] += loc_reaction_forces(2);
901 reactionVec[3] += loc_moment_forces(0);
902 reactionVec[4] += loc_moment_forces(1);
903 reactionVec[5] += loc_moment_forces(2);
911 int nb_integration_pts = data.
getN().size1();
912 auto v = getVolume();
913 auto t_w = getFTensor0IntegrationWeight();
914 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
915 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
916 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
917 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
918 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
919 dataAtPts->wL2DotDotAtPts.clear();
921 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
923 auto piola_scale = dataAtPts->piolaScale;
924 auto alpha_w = alphaW / piola_scale;
925 auto alpha_rho = alphaRho / piola_scale;
927 int nb_base_functions = data.
getN().size2();
931 auto get_ftensor1 = [](
auto &
v) {
936 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
938 auto t_nf = get_ftensor1(nF);
940 for (; bb != nb_dofs / 3; ++bb) {
941 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
942 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
943 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
947 for (; bb != nb_base_functions; ++bb)
961 int nb_integration_pts = getGaussPts().size2();
962 auto v = getVolume();
963 auto t_w = getFTensor0IntegrationWeight();
964 auto t_levi_kirchhoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
965 auto t_omega_grad_dot =
966 getFTensor2FromMat<3, 3>(dataAtPts->rotAxisGradDotAtPts);
967 int nb_base_functions = data.
getN().size2();
973 auto get_ftensor1 = [](
auto &
v) {
977 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
979 auto t_nf = get_ftensor1(nF);
981 for (; bb != nb_dofs / 3; ++bb) {
982 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
984 (
a * alphaOmega) * (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
989 for (; bb != nb_base_functions; ++bb) {
1003 int nb_integration_pts = data.
getN().size1();
1004 auto v = getVolume();
1005 auto t_w = getFTensor0IntegrationWeight();
1007 int nb_base_functions = data.
getN().size2() / 3;
1015 auto get_ftensor1 = [](
auto &
v) {
1022 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1023 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1025 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1027 auto t_nf = get_ftensor1(nF);
1032 for (; bb != nb_dofs / 3; ++bb) {
1033 t_nf(
i) -=
a * (t_row_base_fun(
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2);
1034 t_nf(
i) -=
a * (t_row_base_fun(
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2);
1035 t_nf(
i) +=
a * (t_row_base_fun(
j) *
t_kd(
i,
j));
1040 for (; bb != nb_base_functions; ++bb)
1050 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1052 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1054 auto t_nf = get_ftensor1(nF);
1062 for (; bb != nb_dofs / 3; ++bb) {
1063 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1068 for (; bb != nb_base_functions; ++bb)
1082 int nb_integration_pts = data.
getN().size1();
1083 auto v = getVolume();
1084 auto t_w = getFTensor0IntegrationWeight();
1086 int nb_base_functions = data.
getN().size2() / 9;
1094 auto get_ftensor0 = [](
auto &
v) {
1100 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1101 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1103 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1105 auto t_nf = get_ftensor0(nF);
1110 for (; bb != nb_dofs; ++bb) {
1111 t_nf -=
a * t_row_base_fun(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2;
1112 t_nf -=
a * t_row_base_fun(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2;
1116 for (; bb != nb_base_functions; ++bb) {
1126 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1128 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1130 auto t_nf = get_ftensor0(nF);
1134 t_residuum(
i,
j) = t_h(
i,
j);
1137 for (; bb != nb_dofs; ++bb) {
1138 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1142 for (; bb != nb_base_functions; ++bb) {
1156 int nb_integration_pts = data.
getN().size1();
1157 auto v = getVolume();
1158 auto t_w = getFTensor0IntegrationWeight();
1159 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1160 int nb_base_functions = data.
getN().size2() / 3;
1163 auto get_ftensor1 = [](
auto &
v) {
1168 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1170 auto t_nf = get_ftensor1(nF);
1172 for (; bb != nb_dofs / 3; ++bb) {
1173 double div_row_base = t_row_diff_base_fun(
i,
i);
1174 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1176 ++t_row_diff_base_fun;
1178 for (; bb != nb_base_functions; ++bb) {
1179 ++t_row_diff_base_fun;
1201 for (
auto &ext_strain_block : (*externalStrainVecPtr)) {
1203 if (ext_strain_block.ents.find(fe_ent) != ext_strain_block.ents.end()) {
1206 if (scalingMethodsMap.find(ext_strain_block.blockName) != scalingMethodsMap.end()) {
1207 scale *= scalingMethodsMap.at(ext_strain_block.blockName)->getScale(time);
1210 <<
"No scaling method found for " << ext_strain_block.blockName;
1214 double val =
scale * ext_strain_block.val;
1219 int nb_integration_pts = data.
getN().size1();
1220 auto v = getVolume();
1221 auto t_w = getFTensor0IntegrationWeight();
1223 dataAtPts->externalStrainAtPts.resize(
size_symm, nb_integration_pts,
false);
1224 auto t_external_strain =
1225 getFTensor2SymmetricFromMat<3>(dataAtPts->externalStrainAtPts);
1229 getFTensor4DdgFromPtr<3, 3, 0>(&*dataAtPts->matD.data().begin());
1235 auto get_ftensor2 = [](
auto &
v) {
1237 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1240 int nb_base_functions = data.
getN().size2();
1242 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1243 t_external_strain(
i,
j) = -
t_kd(
i,
j) * val;
1245 auto t_nf = get_ftensor2(nF);
1248 t_T(
i,
j) = t_D(
i,
j,
k,
l) * t_external_strain(
k,
l);
1251 t_residual(
L) =
a * (t_L(
i,
j,
L) * t_T(
i,
j));
1254 for (; bb != nb_dofs / 6; ++bb) {
1255 t_nf(
L) += t_row_base_fun * t_residual(
L);
1259 for (; bb != nb_base_functions; ++bb)
1263 ++t_external_strain;
1275 int nb_integration_pts = getGaussPts().size2();
1278 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1280 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1281 if (tag_length != 9) {
1283 "Number of internal stress components should be 9 but is %d",
1288 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1289 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1290 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1291 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1293 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1294 dataAtPts->internalStressAtPts.clear();
1295 auto t_internal_stress =
1296 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1299 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1300 t_internal_stress(
L) = t_const_stress(
L);
1301 ++t_internal_stress;
1313 int nb_integration_pts = getGaussPts().size2();
1316 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1318 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1319 if (tag_length != 9) {
1321 "Number of internal stress components should be 9 but is %d",
1325 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1328 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1331 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1334 dataAtPts->internalStressAtPts.resize(tag_length, nb_integration_pts,
false);
1335 dataAtPts->internalStressAtPts.clear();
1336 auto t_internal_stress =
1337 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1342 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1343 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1344 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1345 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1349 ++t_internal_stress;
1360 int nb_integration_pts = data.
getN().size1();
1361 auto v = getVolume();
1362 auto t_w = getFTensor0IntegrationWeight();
1367 auto get_ftensor2 = [](
auto &
v) {
1369 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1372 auto t_internal_stress =
1373 getFTensor2FromMat<3, 3>(dataAtPts->internalStressAtPts);
1378 int nb_base_functions = data.
getN().size2();
1380 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1382 auto t_nf = get_ftensor2(nF);
1385 t_symm_stress(
i,
j) =
1386 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1389 t_residual(
L) = t_L(
i,
j,
L) * t_symm_stress(
i,
j);
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;
1411 int nb_integration_pts = data.
getN().size1();
1412 auto v = getVolume();
1413 auto t_w = getFTensor0IntegrationWeight();
1415 auto get_ftensor2 = [](
auto &
v) {
1417 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1420 auto t_internal_stress =
1421 getFTensor1FromMat<9>(dataAtPts->internalStressAtPts);
1428 int nb_base_functions = data.
getN().size2();
1430 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1432 auto t_nf = get_ftensor2(nF);
1435 t_residual(
L) = t_L(M,
L) * t_internal_stress(M);
1438 for (; bb != nb_dofs / 6; ++bb) {
1439 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1443 for (; bb != nb_base_functions; ++bb)
1447 ++t_internal_stress;
1452template <AssemblyType A>
1458 for (
auto &bc : (*bcDispPtr)) {
1460 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1463 int nb_integration_pts = OP::getGaussPts().size2();
1464 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1465 auto t_w = OP::getFTensor0IntegrationWeight();
1466 int nb_base_functions = data.
getN().size2() / 3;
1473 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1475 scale *= scalingMethodsMap.at(bc.blockName)
1478 scale *= scalingMethodsMap.at(bc.blockName)
1479 ->getScale(OP::getFEMethod()->ts_t);
1483 <<
"No scaling method found for " << bc.blockName;
1490 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1491 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1493 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1495 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1499 for (; bb != nb_base_functions; ++bb)
1511 return OP::iNtegrate(data);
1514template <AssemblyType A>
1522 double time = OP::getFEMethod()->ts_t;
1530 for (
auto &bc : (*bcRotPtr)) {
1532 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1534 int nb_integration_pts = OP::getGaussPts().size2();
1535 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1536 auto t_w = OP::getFTensor0IntegrationWeight();
1538 int nb_base_functions = data.
getN().size2() / 3;
1541 auto get_ftensor1 = [](
auto &
v) {
1554 auto get_rotation_angle = [&]() {
1555 double theta = bc.theta;
1556 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1557 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1562 auto get_rotation = [&](
auto theta) {
1564 if (bc.vals.size() == 7) {
1565 t_omega(0) = bc.vals[4];
1566 t_omega(1) = bc.vals[5];
1567 t_omega(2) = bc.vals[6];
1570 t_omega(
i) = OP::getFTensor1Normal()(
i);
1573 t_omega(
i) *= theta;
1580 auto t_R = get_rotation(get_rotation_angle());
1581 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1583 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1585 t_delta(
i) = t_center(
i) - t_coords(
i);
1587 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1589 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1591 for (; bb != nb_dofs / 3; ++bb) {
1592 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1596 for (; bb != nb_base_functions; ++bb)
1609 return OP::iNtegrate(data);
1612template <AssemblyType A>
1616 double time = OP::getFEMethod()->ts_t;
1624 for (
auto &bc : (*bcDispPtr)) {
1626 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1628 auto t_approx_P = getFTensor2FromMat<3, 3>(*piolaStressPtr);
1629 auto t_u = getFTensor1FromMat<3>(*hybridDispPtr);
1630 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1631 auto t_w = OP::getFTensor0IntegrationWeight();
1639 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1640 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1643 <<
"No scaling method found for " << bc.blockName;
1647 double val =
scale * bc.val;
1650 int nb_integration_pts = OP::getGaussPts().size2();
1651 int nb_base_functions = data.
getN().size2();
1653 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1656 t_N(
i) = t_normal(
i);
1660 t_P(
i,
j) = t_N(
i) * t_N(
j);
1665 t_tracton(
i) = t_approx_P(
i,
j) * t_N(
j);
1668 t_res(
i) = t_Q(
i,
j) * t_tracton(
j) + t_P(
i,
j) * 2* t_u(
j) - t_N(
i) * val;
1670 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1672 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1673 t_nf(
i) += (t_w * t_row_base) * t_res(
i);
1677 for (; bb != nb_base_functions; ++bb)
1686 OP::locF *= OP::getMeasure();
1692template <AssemblyType A>
1698 double time = OP::getFEMethod()->ts_t;
1703 int row_nb_dofs = row_data.
getIndices().size();
1704 int col_nb_dofs = col_data.
getIndices().size();
1705 auto &locMat = OP::locMat;
1706 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1712 for (
auto &bc : (*bcDispPtr)) {
1714 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1716 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1717 auto t_w = OP::getFTensor0IntegrationWeight();
1723 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1724 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1727 <<
"No scaling method found for " << bc.blockName;
1730 int nb_integration_pts = OP::getGaussPts().size2();
1731 int row_nb_dofs = row_data.
getIndices().size();
1732 int col_nb_dofs = col_data.
getIndices().size();
1733 int nb_base_functions = row_data.
getN().size2();
1738 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1741 t_N(
i) = t_normal(
i);
1745 t_P(
i,
j) = t_N(
i) * t_N(
j);
1748 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1751 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1752 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1755 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1756 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1763 for (; rr != nb_base_functions; ++rr)
1770 locMat *= OP::getMeasure();
1777template <AssemblyType A>
1783 double time = OP::getFEMethod()->ts_t;
1788 int row_nb_dofs = row_data.
getIndices().size();
1789 int col_nb_dofs = col_data.
getIndices().size();
1790 auto &locMat = OP::locMat;
1791 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1797 for (
auto &bc : (*bcDispPtr)) {
1799 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1801 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1802 auto t_w = OP::getFTensor0IntegrationWeight();
1811 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1812 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1815 <<
"No scaling method found for " << bc.blockName;
1818 int nb_integration_pts = OP::getGaussPts().size2();
1819 int nb_base_functions = row_data.
getN().size2();
1822 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1825 t_N(
i) = t_normal(
i);
1829 t_P(
i,
j) = t_N(
i) * t_N(
j);
1834 t_d_res(
i,
j) = t_Q(
i,
j);
1837 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1838 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1841 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1843 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
1850 for (; rr != nb_base_functions; ++rr)
1857 locMat *= OP::getMeasure();
1864 return OP::iNtegrate(data);
1869 return OP::iNtegrate(row_data, col_data);
1874 return OP::iNtegrate(row_data, col_data);
1877template <AssemblyType A>
1881 double time = OP::getFEMethod()->ts_t;
1889 for (
auto &bc : (*bcDispPtr)) {
1891 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1896 auto [block_name, v_analytical_expr] =
1901 int nb_integration_pts = OP::getGaussPts().size2();
1902 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1903 auto t_w = OP::getFTensor0IntegrationWeight();
1904 int nb_base_functions = data.
getN().size2() / 3;
1906 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
1912 auto t_bc_disp = getFTensor1FromMat<3>(v_analytical_expr);
1914 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1915 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1918 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1920 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1924 for (; bb != nb_base_functions; ++bb)
1938 return OP::iNtegrate(data);
1947 int nb_integration_pts = getGaussPts().size2();
1948 int nb_base_functions = data.
getN().size2();
1950 double time = getFEMethod()->ts_t;
1956 if (this->locF.size() != nb_dofs)
1958 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
1961 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
1964 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
1966 auto t_w = getFTensor0IntegrationWeight();
1967 auto t_coords = getFTensor1CoordsAtGaussPts();
1969 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
1971 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1973 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
1974 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
1976 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1977 t_f(
i) -= time_scale * t_w * t_row_base * tau * (t_val(
i) *
scale);
1982 for (; rr != nb_base_functions; ++rr)
1989 this->locF *= 2. * getMeasure();
1995 for (
auto &bc : *(bcData)) {
1996 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1998 double time_scale = 1;
1999 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2000 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2006 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
2010 return -y * (y - 1) / 0.25;
2012 CHKERR integrate_rhs(bc, calc_tau, time_scale);
2015 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2029 int nb_integration_pts = getGaussPts().size2();
2030 int nb_base_functions = data.
getN().size2();
2032 double time = getFEMethod()->ts_t;
2038 if (this->locF.size() != nb_dofs)
2040 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2043 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2048 auto t_w = getFTensor0IntegrationWeight();
2049 auto t_coords = getFTensor1CoordsAtGaussPts();
2050 auto t_normal = getFTensor1NormalsAtGaussPts();
2052 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2054 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2055 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2058 (time_scale * t_w * tau *
scale * val) * t_normal(
i) / t_normal.l2();
2060 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2062 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2063 t_f(
i) -= t_row_base * t_val(
i);
2068 for (; rr != nb_base_functions; ++rr)
2076 this->locF *= 2. * getMeasure();
2082 for (
auto &bc : *(bcData)) {
2083 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2085 double time_scale = 1;
2086 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2087 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2093 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2106 int nb_integration_pts = getGaussPts().size2();
2107 int nb_base_functions = data.
getN().size2();
2109 double time = getFEMethod()->ts_t;
2115 if (this->locF.size() != nb_dofs)
2117 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2122 for (
auto &bc : *(bcData)) {
2123 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2127 auto [block_name, v_analytical_expr] =
2130 getFTensor1FromMat<3>(v_analytical_expr);
2132 auto t_w = getFTensor0IntegrationWeight();
2133 auto t_coords = getFTensor1CoordsAtGaussPts();
2135 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2137 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2139 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2141 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2142 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2147 for (; rr != nb_base_functions; ++rr)
2155 this->locF *= 2. * getMeasure();
2164 int nb_integration_pts = row_data.
getN().size1();
2165 int row_nb_dofs = row_data.
getIndices().size();
2166 int col_nb_dofs = col_data.
getIndices().size();
2167 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2169 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2172 auto v = getVolume();
2173 auto t_w = getFTensor0IntegrationWeight();
2174 int row_nb_base_functions = row_data.
getN().size2();
2176 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2179 for (; rr != row_nb_dofs / 3; ++rr) {
2181 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2182 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2183 double div_col_base = t_col_diff_base_fun(
i,
i);
2184 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2186 ++t_col_diff_base_fun;
2190 for (; rr != row_nb_base_functions; ++rr)
2201 if (alphaW < std::numeric_limits<double>::epsilon() &&
2202 alphaRho < std::numeric_limits<double>::epsilon())
2205 const int nb_integration_pts = row_data.
getN().size1();
2206 const int row_nb_dofs = row_data.
getIndices().size();
2207 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2209 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2215 auto v = getVolume();
2216 auto t_w = getFTensor0IntegrationWeight();
2218 auto piola_scale = dataAtPts->piolaScale;
2219 auto alpha_w = alphaW / piola_scale;
2220 auto alpha_rho = alphaRho / piola_scale;
2222 int row_nb_base_functions = row_data.
getN().size2();
2225 double ts_scale = alpha_w * getTSa();
2226 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2227 ts_scale += alpha_rho * getTSaa();
2229 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2230 double a =
v * t_w * ts_scale;
2233 for (; rr != row_nb_dofs / 3; ++rr) {
2236 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2237 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2238 const double b =
a * t_row_base_fun * t_col_base_fun;
2247 for (; rr != row_nb_base_functions; ++rr)
2266 int nb_integration_pts = row_data.
getN().size1();
2267 int row_nb_dofs = row_data.
getIndices().size();
2268 int col_nb_dofs = col_data.
getIndices().size();
2269 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2272 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2274 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2276 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2278 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2280 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2282 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2285 auto v = getVolume();
2286 auto t_w = getFTensor0IntegrationWeight();
2288 int row_nb_base_functions = row_data.
getN().size2();
2291 auto t_approx_P_adjoint_log_du_dP =
2292 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2294 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2297 for (; rr != row_nb_dofs / 6; ++rr) {
2300 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2302 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2304 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2312 for (; rr != row_nb_base_functions; ++rr)
2315 ++t_approx_P_adjoint_log_du_dP;
2331 int nb_integration_pts = row_data.
getN().size1();
2332 int row_nb_dofs = row_data.
getIndices().size();
2333 int col_nb_dofs = col_data.
getIndices().size();
2334 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2336 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2340 auto v = getVolume();
2341 auto t_w = getFTensor0IntegrationWeight();
2344 int row_nb_base_functions = row_data.
getN().size2();
2346 auto t_approx_P_adjoint_log_du_dP =
2347 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
2349 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2352 for (; rr != row_nb_dofs / 6; ++rr) {
2353 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2354 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2355 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2357 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2364 for (; rr != row_nb_base_functions; ++rr)
2367 ++t_approx_P_adjoint_log_du_dP;
2379 int row_nb_dofs = row_data.
getIndices().size();
2380 int col_nb_dofs = col_data.
getIndices().size();
2381 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2384 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2386 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2388 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2390 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2392 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2394 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2404 auto v = getVolume();
2405 auto t_w = getFTensor0IntegrationWeight();
2406 auto t_approx_P_adjoint_log_du_domega =
2407 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2409 int row_nb_base_functions = row_data.
getN().size2();
2412 int nb_integration_pts = row_data.
getN().size1();
2414 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2418 for (; rr != row_nb_dofs / 6; ++rr) {
2420 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2421 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2422 double v =
a * t_row_base_fun * t_col_base_fun;
2423 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2430 for (; rr != row_nb_base_functions; ++rr)
2434 ++t_approx_P_adjoint_log_du_domega;
2443 int nb_integration_pts = getGaussPts().size2();
2444 int row_nb_dofs = row_data.
getIndices().size();
2445 int col_nb_dofs = col_data.
getIndices().size();
2446 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2450 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2451 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2453 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2454 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2456 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2457 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2465 auto v = getVolume();
2466 auto t_w = getFTensor0IntegrationWeight();
2467 auto t_levi_kirchhoff_du = getFTensor2FromMat<3, size_symm>(
2468 dataAtPts->leviKirchhoffdLogStreatchAtPts);
2469 int row_nb_base_functions = row_data.
getN().size2();
2471 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2474 for (; rr != row_nb_dofs / 3; ++rr) {
2475 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2476 const double b =
a * t_row_base_fun;
2478 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2479 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2485 for (; rr != row_nb_base_functions; ++rr) {
2489 ++t_levi_kirchhoff_du;
2505 int nb_integration_pts = getGaussPts().size2();
2506 int row_nb_dofs = row_data.
getIndices().size();
2507 int col_nb_dofs = col_data.
getIndices().size();
2508 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2511 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2513 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2515 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2520 auto v = getVolume();
2521 auto t_w = getFTensor0IntegrationWeight();
2523 int row_nb_base_functions = row_data.
getN().size2();
2525 auto t_levi_kirchhoff_dP =
2526 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2528 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2531 for (; rr != row_nb_dofs / 3; ++rr) {
2532 double b =
a * t_row_base_fun;
2534 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2535 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2536 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2542 for (; rr != row_nb_base_functions; ++rr) {
2547 ++t_levi_kirchhoff_dP;
2555 int nb_integration_pts = getGaussPts().size2();
2556 int row_nb_dofs = row_data.
getIndices().size();
2557 int col_nb_dofs = col_data.
getIndices().size();
2559 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2561 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2568 auto v = getVolume();
2569 auto t_w = getFTensor0IntegrationWeight();
2570 auto t_levi_kirchoff_dP =
2571 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2573 int row_nb_base_functions = row_data.
getN().size2();
2576 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2579 for (; rr != row_nb_dofs / 3; ++rr) {
2580 double b =
a * t_row_base_fun;
2581 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2582 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2583 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2584 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2591 for (; rr != row_nb_base_functions; ++rr) {
2595 ++t_levi_kirchoff_dP;
2603 int nb_integration_pts = getGaussPts().size2();
2604 int row_nb_dofs = row_data.
getIndices().size();
2605 int col_nb_dofs = col_data.
getIndices().size();
2606 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2609 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2611 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2613 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2626 auto v = getVolume();
2627 auto ts_a = getTSa();
2628 auto t_w = getFTensor0IntegrationWeight();
2629 auto t_levi_kirchhoff_domega =
2630 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffdOmegaAtPts);
2631 int row_nb_base_functions = row_data.
getN().size2();
2634 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2636 double c =
a * alphaOmega * ts_a;
2638 for (; rr != row_nb_dofs / 3; ++rr) {
2639 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2640 const double b =
a * t_row_base_fun;
2643 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2644 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2645 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2653 for (; rr != row_nb_base_functions; ++rr) {
2658 ++t_levi_kirchhoff_domega;
2666 if (dataAtPts->matInvD.size1() ==
size_symm &&
2667 dataAtPts->matInvD.size2() ==
size_symm) {
2680 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2683 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2685 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2687 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
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();
2696 auto v = getVolume();
2697 auto t_w = getFTensor0IntegrationWeight();
2698 int row_nb_base_functions = row_data.
getN().size2() / 3;
2705 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2706 &*dataAtPts->matInvD.data().begin());
2709 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2713 for (; rr != row_nb_dofs / 3; ++rr) {
2715 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2716 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2717 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2725 for (; rr != row_nb_base_functions; ++rr)
2738 if (dataAtPts->matInvD.size1() ==
size_symm &&
2739 dataAtPts->matInvD.size2() ==
size_symm) {
2753 int nb_integration_pts = getGaussPts().size2();
2754 int row_nb_dofs = row_data.
getIndices().size();
2755 int col_nb_dofs = col_data.
getIndices().size();
2757 auto v = getVolume();
2758 auto t_w = getFTensor0IntegrationWeight();
2759 int row_nb_base_functions = row_data.
getN().size2() / 9;
2766 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2767 &*dataAtPts->matInvD.data().begin());
2770 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2774 for (; rr != row_nb_dofs; ++rr) {
2776 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2778 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2785 for (; rr != row_nb_base_functions; ++rr)
2796 if (dataAtPts->matInvD.size1() ==
size_symm &&
2797 dataAtPts->matInvD.size2() ==
size_symm) {
2811 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2814 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
2819 int nb_integration_pts = getGaussPts().size2();
2820 int row_nb_dofs = row_data.
getIndices().size();
2821 int col_nb_dofs = col_data.
getIndices().size();
2823 auto v = getVolume();
2824 auto t_w = getFTensor0IntegrationWeight();
2825 int row_nb_base_functions = row_data.
getN().size2() / 9;
2834 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
2835 &*dataAtPts->matInvD.data().begin());
2838 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2841 auto t_m = get_ftensor1(
K, 0, 0);
2844 for (; rr != row_nb_dofs; ++rr) {
2846 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2847 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2855 for (; rr != row_nb_base_functions; ++rr)
2874 int nb_integration_pts = row_data.
getN().size1();
2875 int row_nb_dofs = row_data.
getIndices().size();
2876 int col_nb_dofs = col_data.
getIndices().size();
2878 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2881 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2883 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2885 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2890 auto v = getVolume();
2891 auto t_w = getFTensor0IntegrationWeight();
2892 int row_nb_base_functions = row_data.
getN().size2() / 3;
2895 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2896 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2900 for (; rr != row_nb_dofs / 3; ++rr) {
2903 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2906 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2907 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2908 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
2916 for (; rr != row_nb_base_functions; ++rr)
2936 int nb_integration_pts = row_data.
getN().size1();
2937 int row_nb_dofs = row_data.
getIndices().size();
2938 int col_nb_dofs = col_data.
getIndices().size();
2940 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2942 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
2945 auto v = getVolume();
2946 auto t_w = getFTensor0IntegrationWeight();
2947 int row_nb_base_functions = row_data.
getN().size2() / 9;
2950 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2951 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2955 for (; rr != row_nb_dofs; ++rr) {
2958 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
2961 auto t_m = get_ftensor2(
K, rr, 0);
2962 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2963 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
2971 for (; rr != row_nb_base_functions; ++rr)
2984 if (tagSense != getSkeletonSense())
2987 auto get_tag = [&](
auto name) {
2988 auto &mob = getPtrFE()->mField.get_moab();
2994 auto get_tag_value = [&](
auto &&tag,
int dim) {
2995 auto &mob = getPtrFE()->mField.get_moab();
2996 auto face = getSidePtrFE()->getFEEntityHandle();
2997 std::vector<double> value(dim);
2998 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3002 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3003 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3005 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3006 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3011 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3012 Tag th_detF = create_tag(
"detF", 1);
3013 Tag th_traction = create_tag(
"traction", 3);
3014 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3016 Tag th_energy = create_tag(
"Energy", 1);
3018 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
3019 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
3020 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
3022 auto t_normal = getFTensor1NormalsAtGaussPts();
3023 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
3033 if (dataAtPts->energyAtPts.size() == 0) {
3035 dataAtPts->energyAtPts.resize(getGaussPts().size2());
3036 dataAtPts->energyAtPts.clear();
3045 auto set_float_precision = [](
const double x) {
3046 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3053 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3055 v = set_float_precision(
v);
3056 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3063 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3066 for (
auto &
a :
v.data())
3067 a = set_float_precision(
a);
3068 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3069 &*
v.data().begin());
3077 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3079 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3081 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3083 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3085 t_m(
i,
j) = t_d(
i,
j);
3086 for (
auto &
v :
m.data())
3087 v = set_float_precision(
v);
3088 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3089 &*
m.data().begin());
3093 const auto nb_gauss_pts = getGaussPts().size2();
3094 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3097 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3099 CHKERR save_vec_tag(th_traction, t_traction, gg);
3101 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3102 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3103 CHKERR save_scal_tag(th_energy, t_energy, gg);
3107 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3108 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3109 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3118 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3119 std::vector<FieldSpace> spaces, std::string geom_field_name,
3120 boost::shared_ptr<Range> crack_front_edges_ptr) {
3123 constexpr bool scale_l2 =
false;
3133 boost::shared_ptr<Range> edges_ptr)
3142 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3145 return OP::doWork(side, type, data);
3150 boost::shared_ptr<Range> edgesPtr;
3153 auto jac = boost::make_shared<MatrixDouble>();
3154 auto det = boost::make_shared<VectorDouble>();
3157 ? crack_front_edges_ptr
3158 : boost::make_shared<
Range>()));
3171 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3172 std::vector<FieldSpace> spaces, std::string geom_field_name,
3173 boost::shared_ptr<Range> crack_front_edges_ptr) {
3176 constexpr bool scale_l2 =
false;
3180 "Scale L2 Ainsworth Legendre base is not implemented");
3189 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3190 std::vector<FieldSpace> spaces, std::string geom_field_name,
3191 boost::shared_ptr<Range> crack_front_edges_ptr) {
3195 auto jac = boost::make_shared<MatrixDouble>();
3196 auto det = boost::make_shared<VectorDouble>();
3198 geom_field_name, jac));
3204 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3206 if (scale_l2_ainsworth_legendre_base) {
3214 boost::shared_ptr<MatrixDouble> jac,
3215 boost::shared_ptr<Range> edges_ptr)
3224 if (type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3227 return OP::doWork(side, type, data);
3232 boost::shared_ptr<Range> edgesPtr;
3235 auto jac = boost::make_shared<MatrixDouble>();
3236 auto det = boost::make_shared<VectorDouble>();
3238 geom_field_name, jac,
3240 : boost::make_shared<
Range>()));
3247 nullptr,
nullptr,
nullptr);
3262 dataAtPts->faceMaterialForceAtPts.resize(3, getGaussPts().size2(),
false);
3263 dataAtPts->normalPressureAtPts.resize(getGaussPts().size2(),
false);
3264 if (getNinTheLoop() == 0) {
3265 dataAtPts->faceMaterialForceAtPts.clear();
3266 dataAtPts->normalPressureAtPts.clear();
3268 auto loop_size = getLoopSize();
3269 if (loop_size == 1) {
3270 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3271 auto pstatus = numebered_fe_ptr->getPStatus();
3272 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3279 auto t_normal = getFTensor1NormalsAtGaussPts();
3280 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3281 dataAtPts->faceMaterialForceAtPts);
3284 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
3285 auto t_grad_u_gamma =
3286 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->gradHybridDispAtPts);
3288 getFTensor2SymmetricFromMat<SPACE_DIM>(dataAtPts->logStretchTensorAtPts);
3292 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3294 t_N(
I) = t_normal(
I);
3299 (t_grad_u_gamma(
i,
j) - t_grad_u_gamma(
j,
i)) / 2. + t_strain(
i,
j);
3300 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3302 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2) / loop_size;
3305 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3317 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3319 t_N(
I) = t_normal(
I);
3322 t_T(
I) += t_N(
J) * (t_grad_u_gamma(
i,
I) * t_P(
i,
J)) / loop_size;
3325 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3338 "Grffith energy release "
3339 "selector not implemented");
3343 auto side_fe_ptr = getSidePtrFE();
3344 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3345 auto pstatus = side_fe_mi_ptr->getPStatus();
3347 auto owner = side_fe_mi_ptr->getOwnerProc();
3349 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3350 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3351 << getNinTheLoop() <<
" loop size " << getLoopSize();
3363 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3364 auto pstatus = fe_mi_ptr->getPStatus();
3366 auto owner = fe_mi_ptr->getOwnerProc();
3368 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3369 <<
" " << getPtrFE()->mField.get_comm_rank();
3377 double face_pressure = 0.;
3378 auto t_T = getFTensor1FromMat<SPACE_DIM>(
3379 dataAtPts->faceMaterialForceAtPts);
3382 auto t_w = getFTensor0IntegrationWeight();
3383 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3384 t_face_T(
I) += t_w * t_T(
I);
3385 face_pressure += t_w * t_p;
3390 t_face_T(
I) *= getMeasure();
3391 face_pressure *= getMeasure();
3393 auto get_tag = [&](
auto name,
auto dim) {
3394 auto &moab = getPtrFE()->mField.get_moab();
3396 double def_val[] = {0., 0., 0.};
3397 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3398 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3403 auto set_tag = [&](
auto &&tag,
auto ptr) {
3404 auto &moab = getPtrFE()->mField.get_moab();
3405 auto face = getPtrFE()->getFEEntityHandle();
3406 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3409 set_tag(get_tag(
"MaterialForce", 3), &t_face_T(0));
3410 set_tag(get_tag(
"FacePressure", 1), &face_pressure);
3416#ifdef ENABLE_PYTHON_BINDING
3418 MoFEMErrorCode AnalyticalExprPython::analyticalExprInit(
const std::string py_file) {
3423 auto main_module = bp::import(
"__main__");
3424 mainNamespace = main_module.attr(
"__dict__");
3425 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
3427 analyticalDispFun = mainNamespace[
"analytical_disp"];
3428 analyticalTractionFun = mainNamespace[
"analytical_traction"];
3430 }
catch (bp::error_already_set
const &) {
3440 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3441 const std::string& block_name, np::ndarray &analytical_expr
3448 analytical_expr = bp::extract<np::ndarray>(
3449 analyticalDispFun(delta_t,
t, x, y, z, block_name));
3451 }
catch (bp::error_already_set
const &) {
3461 double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
3462 const std::string& block_name, np::ndarray &analytical_expr
3469 analytical_expr = bp::extract<np::ndarray>(
3470 analyticalTractionFun(delta_t,
t, x, y, z, block_name));
3472 }
catch (bp::error_already_set
const &) {
3480boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
3482inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
3484 auto dtype = np::dtype::get_builtin<double>();
3485 auto size = bp::make_tuple(nb_gauss_pts);
3486 auto stride = bp::make_tuple(3 *
sizeof(
double));
3487 return (np::from_data(&data[
id], dtype, size, stride, bp::object()));
3496 const std::string block_name) {
3498#ifdef ENABLE_PYTHON_BINDING
3499 if (
auto analytical_expr_ptr = AnalyticalExprPythonWeakPtr.lock()) {
3503 bp::list python_coords;
3505 for (
int idx = 0; idx < 3; ++idx) {
3506 python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
3509 np::ndarray np_analytical_expr = np::empty(bp::make_tuple(nb_gauss_pts, 3),
3510 np::dtype::get_builtin<double>());
3512 auto disp_block_name =
"(.*)ANALYTICAL_DISPLACEMENT(.*)";
3513 auto traction_block_name =
"(.*)ANALYTICAL_TRACTION(.*)";
3515 std::regex reg_disp_name(disp_block_name);
3516 std::regex reg_traction_name(traction_block_name);
3517 if (std::regex_match(block_name, reg_disp_name)) {
3519 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3520 bp::extract<np::ndarray>(python_coords[1]),
3521 bp::extract<np::ndarray>(python_coords[2]), block_name,
3522 np_analytical_expr),
3523 "Failed analytical_disp() python call");
3525 else if (std::regex_match(block_name, reg_traction_name)) {
3527 delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
3528 bp::extract<np::ndarray>(python_coords[1]),
3529 bp::extract<np::ndarray>(python_coords[2]), block_name,
3530 np_analytical_expr),
3531 "Failed analytical_traction() python call");
3534 "Unknown analytical block");
3537 double *analytical_expr_val_ptr =
3538 reinterpret_cast<double *
>(np_analytical_expr.get_data());
3541 v_analytical_expr.resize(3, nb_gauss_pts,
false);
3542 for (
size_t gg = 0; gg < nb_gauss_pts; ++gg) {
3543 for (
int idx = 0; idx < 3; ++idx)
3544 v_analytical_expr(idx, gg) = *(analytical_expr_val_ptr + (3 * gg + idx));
3546 return v_analytical_expr;
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.
auto sort_eigen_vals(A &eig, B &eigen_vec)
MatrixDouble analytical_expr_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_ref_coords, const std::string block_name)
[Analytical displacement function from python]
static constexpr auto size_symm
auto getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
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
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 dofs values
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 dofs on entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Scale base functions by inverses of measure of element.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &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)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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)
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)
