14#include <boost/math/constants/constants.hpp>
32 int nb_integration_pts = getGaussPts().size2();
34 auto t_P = dataAtPts->getFTensorApproxP(getGaussPts().size2());
35 auto t_F = dataAtPts->getFTensorSmallH(getGaussPts().size2());
38 auto get_eshelby_stress =
40 DL>::size(dataAtPts->SigmaAtPts, nb_integration_pts);
41 auto t_eshelby_stress = get_eshelby_stress();
45 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
46 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
62 int nb_integration_pts = getGaussPts().size2();
79 *dataAtPts->getStretchTensorAtPts(), nb_integration_pts);
81 *dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
83 *dataAtPts->getStretchH1AtPts(), nb_integration_pts);
85 *dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
87 *dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
89 *dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
91 *dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
93 *dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
96 *dataAtPts->getDeformationGradient(), nb_integration_pts);
98 dataAtPts->hdOmegaAtPts, nb_integration_pts);
100 dataAtPts->hdLogStretchAtPts, nb_integration_pts);
103 dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
105 dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
107 dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
109 dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
112 dataAtPts->rotMatAtPts, nb_integration_pts);
114 *dataAtPts->getEigenVals(), nb_integration_pts);
116 *dataAtPts->getEigenVecs(), nb_integration_pts);
117 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
119 dataAtPts->eigenValsC, nb_integration_pts);
121 dataAtPts->eigenVecsC, nb_integration_pts);
122 dataAtPts->nbUniqC.resize(nb_integration_pts,
false);
125 dataAtPts->logStretch2H1AtPts, nb_integration_pts);
127 dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
130 dataAtPts->internalStressAtPts, nb_integration_pts);
131 dataAtPts->internalStressAtPts.clear();
134 auto t_h = dataAtPts->getFTensorSmallH(getGaussPts().size2());
135 auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(getGaussPts().size2());
137 dataAtPts->getFTensorSmallHdLogStretch(getGaussPts().size2());
138 auto t_levi_kirchhoff =
139 dataAtPts->getFTensorLeviKirchhoff(getGaussPts().size2());
140 auto t_levi_kirchhoff_domega =
141 dataAtPts->getFTensorLeviKirchhoffdOmega(getGaussPts().size2());
142 auto t_levi_kirchhoff_dstreach =
143 dataAtPts->getFTensorLeviKirchhoffdLogStretch(getGaussPts().size2());
144 auto t_levi_kirchhoff_dP =
145 dataAtPts->getFTensorLeviKirchhoffP(getGaussPts().size2());
146 auto t_approx_P_adjoint_dstretch =
147 dataAtPts->getFTensorAdjointPdstretch(getGaussPts().size2());
148 auto t_approx_P_adjoint_log_du =
149 dataAtPts->getFTensorAdjointPdU(getGaussPts().size2());
150 auto t_approx_P_adjoint_log_du_dP =
151 dataAtPts->getFTensorAdjointPdUdP(getGaussPts().size2());
152 auto t_approx_P_adjoint_log_du_domega =
153 dataAtPts->getFTensorAdjointPdUdOmega(getGaussPts().size2());
154 auto t_R = dataAtPts->getFTensorRotMat(getGaussPts().size2());
155 auto t_u = dataAtPts->getFTensorStretch(getGaussPts().size2());
156 auto t_diff_u = dataAtPts->getFTensorDiffStretch(getGaussPts().size2());
157 auto t_eigen_vals = dataAtPts->getFTensorEigenVals(getGaussPts().size2());
158 auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(getGaussPts().size2());
159 auto &nbUniq = dataAtPts->nbUniq;
162 auto t_eigen_vals_C = dataAtPts->getFTensorEigenValsC(nb_integration_pts);
163 auto t_eigen_vecs_C = dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
164 auto &nbUniqC = dataAtPts->nbUniqC;
168 auto t_u_h1 = dataAtPts->getFTensorStretchH1(getGaussPts().size2());
169 auto t_diff_u_h1 = dataAtPts->getFTensorDiffStretchH1(getGaussPts().size2());
170 auto t_log_stretch_total =
171 dataAtPts->getFTensorLogStretchTotal(getGaussPts().size2());
172 auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(getGaussPts().size2());
175 auto t_grad_h1 = dataAtPts->getFTensorSmallWGradH1(getGaussPts().size2());
176 auto t_omega = dataAtPts->getFTensorRotAxis(getGaussPts().size2());
177 auto t_approx_P = dataAtPts->getFTensorApproxP(getGaussPts().size2());
178 auto t_log_u = dataAtPts->getFTensorLogStretch(getGaussPts().size2());
181 auto t_omega0 = dataAtPts->getFTensorRotAxis0(getGaussPts().size2());
189 ++t_levi_kirchhoff_domega;
190 ++t_levi_kirchhoff_dstreach;
191 ++t_levi_kirchhoff_dP;
192 ++t_approx_P_adjoint_dstretch;
193 ++t_approx_P_adjoint_log_du;
194 ++t_approx_P_adjoint_log_du_dP;
195 ++t_approx_P_adjoint_log_du_domega;
208 ++t_log_stretch_total;
220 auto bound_eig = [&](
auto &eig) {
222 const auto zero = std::exp(std::numeric_limits<double>::min_exponent);
223 const auto large = std::exp(std::numeric_limits<double>::max_exponent);
224 for (
int ii = 0; ii != 3; ++ii) {
225 eig(ii) = std::min(std::max(zero, eig(ii)), large);
230 auto calculate_log_stretch = [&]() {
234 eigen_vec(
i,
j) = t_log_u(
i,
j);
236 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
240 t_nb_uniq = get_uniq_nb<3>(&eig(0));
244 t_eigen_vals(
i) = eig(
i);
245 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
248 auto get_t_diff_u = [&]() {
253 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
255 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
260 auto calculate_total_stretch = [&](
auto &t_h1) {
264 t_log_u2_h1(
i,
j) = 0;
265 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
273 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
274 t_eigen_vec(
i,
j) = t_C_h1(
i,
j);
276 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
279 t_nb_uniq_C = get_uniq_nb<3>(&t_eig(0));
280 if (t_nb_uniq_C < 3) {
286 t_eigen_vals_C(
i) = t_eig(
i);
287 t_eigen_vecs_C(
i,
j) = t_eigen_vec(
i,
j);
291 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) / 2 + t_log_u(
i,
j);
296 auto no_h1_loop = [&]() {
305 "no_h1_loop is implemented only for LARGE_ROT");
308 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
316 CHKERR calculate_log_stretch();
318 CHKERR calculate_total_stretch(t_h1);
320 t_u_h1(
i,
j) = t_u(
i,
j);
321 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i,
j,
k,
l);
323 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
328 auto large_rot = [&]() {
332 t_diff_diff_R(
i,
j,
k,
l) =
335 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
337 t_approx_P_adjoint_dstretch(
l,
k) = t_R(
i,
l) * t_approx_P(
i,
k);
338 t_approx_P_adjoint_log_du(
L) =
339 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
341 t_levi_kirchhoff(
m) =
342 t_diff_R(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k));
345 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
346 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
348 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
349 t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
352 t_A(
k,
l,
m) = t_diff_R(
i,
l,
m) * t_approx_P(
i,
k);
353 t_approx_P_adjoint_log_du_domega(
m,
L) =
354 t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L);
356 t_levi_kirchhoff_dstreach(
m,
L) =
357 t_diff_R(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k));
358 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
359 t_levi_kirchhoff_domega(
m,
n) =
360 t_diff_diff_R(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k));
364 auto moderate_rot = [&]() {
367 t_delta_omega(
m) = t_omega(
m) - t_omega0(
m);
371 t_R(
i,
j) = t_R0(
i,
j);
374 (t_u(
l,
k) + levi_civita(
l,
k,
m) * t_delta_omega(
m));
376 t_approx_P_adjoint_dstretch(
l,
k) = t_R0(
i,
l) * t_approx_P(
i,
k);
377 t_approx_P_adjoint_log_du(
L) =
378 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
380 t_levi_kirchhoff(
m) = t_R0(
i,
l) * levi_civita(
l,
k,
m) *
384 t_h_domega(
i,
k,
m) = t_R0(
i,
l) * levi_civita(
l,
k,
m);
385 t_h_dlog_u(
i,
k,
L) = t_R0(
i,
l) * t_Ldiff_u(
l,
k,
L);
387 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
388 t_R0(
i,
l) * t_Ldiff_u(
l,
k,
L);
389 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
391 t_levi_kirchhoff_dstreach(
m,
L) = 0;
392 t_levi_kirchhoff_dP(
m,
i,
k) =
393 t_R0(
i,
l) * levi_civita(
l,
k,
m);
394 t_levi_kirchhoff_domega(
m,
n) = 0;
408 "rotationSelector not handled");
417 auto large_loop = [&]() {
427 "rotSelector should be large or small");
430 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
441 "Selected grad approximator not handled");
445 CHKERR calculate_log_stretch();
447 CHKERR calculate_total_stretch(t_h1);
449 t_u_h1(
l,
k) = t_u(
l, o) * t_h1(o,
k);
450 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i, o,
k,
l) * t_h1(o,
j);
452 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
460 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
461 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
462 t_diff_diff_R(
i,
j,
l,
m) = 0;
468 t_diff_diff_R(
i,
j,
k,
l) =
474 "rotationSelector not handled");
478 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
481 t_approx_P_adjoint_dstretch(
l, o) =
482 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
483 t_approx_P_adjoint_log_du(
L) =
484 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
487 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
491 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
492 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
494 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
495 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
498 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
499 t_approx_P_adjoint_log_du_domega(
m,
L) =
500 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
502 t_levi_kirchhoff_dstreach(
m,
L) =
503 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
505 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
506 t_levi_kirchhoff_domega(
m,
n) =
507 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
516 auto moderate_loop = [&]() {
526 "rotSelector should be large or small");
529 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
540 "Selected grad approximator not handled");
544 CHKERR calculate_log_stretch();
546 CHKERR calculate_total_stretch(t_h1);
550 t_u_h1(
l,
k) = (
t_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
551 t_diff_u_h1(
i,
j,
k,
l) = t_diff(
i, o,
k,
l) * t_h1(o,
j);
553 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
561 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
562 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
563 t_diff_diff_R(
i,
j,
l,
m) = 0;
569 t_diff_diff_R(
i,
j,
k,
l) =
575 "rotationSelector not handled");
579 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
582 t_approx_P_adjoint_dstretch(
l, o) =
583 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
584 t_approx_P_adjoint_log_du(
L) =
585 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
588 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
592 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
593 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
595 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
596 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
599 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
600 t_approx_P_adjoint_log_du_domega(
m,
L) =
601 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
603 t_levi_kirchhoff_dstreach(
m,
L) =
604 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
606 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
607 t_levi_kirchhoff_domega(
m,
n) =
608 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
617 auto small_loop = [&]() {
624 "rotSelector should be small");
627 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
636 "gradApproximator not handled");
642 "stretchSelector should be linear for small loop");
644 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
645 t_u_h1(
i,
j) = t_u(
i,
j);
646 t_diff_u_h1(
i,
j,
k,
l) =
648 t_diff_u_h1(
i,
j,
k,
l) /= 2.;
649 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
651 t_log_u2_h1(
i,
j) = 0;
652 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
654 t_R(
i,
j) =
t_kd(
i,
j) + levi_civita(
i,
j,
k) * t_omega(
k);
655 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
657 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
658 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
661 t_approx_P_adjoint_dstretch(
i,
j) = t_approx_P(
i,
j);
662 t_approx_P_adjoint_log_du(
L) =
663 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
664 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
665 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
668 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
669 t_levi_kirchhoff_dstreach(
m,
L) = 0;
670 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
671 t_levi_kirchhoff_domega(
m,
n) = 0;
697 "gradApproximator not handled");
710 auto N_InLoop = getNinTheLoop();
711 auto sensee = getSkeletonSense();
712 auto nb_gauss_pts = getGaussPts().size2();
713 auto t_normal = getFTensor1NormalsAtGaussPts();
715 auto t_sigma_ptr = dataAtPts->getFTensorApproxP(getGaussPts().size2());
718 dataAtPts->tractionAtPts, nb_gauss_pts);
720 dataAtPts->tractionAtPts.clear();
723 auto t_traction = get_tracion();
724 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
725 t_traction(
i) = t_sigma_ptr(
i,
j) * sensee * (t_normal(
j) / t_normal.l2());
737 if (blockEntities.find(getFEEntityHandle()) == blockEntities.end()) {
743 int nb_integration_pts = getGaussPts().size2();
744 auto t_w = getFTensor0IntegrationWeight();
745 auto t_traction = dataAtPts->getFTensorTraction(nb_integration_pts);
746 auto t_coords = getFTensor1CoordsAtGaussPts();
747 auto t_spatial_disp = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
755 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
756 loc_reaction_forces(
i) += (t_traction(
i)) * t_w * getMeasure();
757 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
759 loc_moment_forces(
i) +=
760 (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j)) *
761 t_traction(
k) * t_w * getMeasure();
768 reactionVec[0] += loc_reaction_forces(0);
769 reactionVec[1] += loc_reaction_forces(1);
770 reactionVec[2] += loc_reaction_forces(2);
771 reactionVec[3] += loc_moment_forces(0);
772 reactionVec[4] += loc_moment_forces(1);
773 reactionVec[5] += loc_moment_forces(2);
781 int nb_integration_pts = data.
getN().size1();
782 auto v = getVolume();
783 auto t_w = getFTensor0IntegrationWeight();
784 auto t_div_P = dataAtPts->getFTensorDivP(nb_integration_pts);
785 auto t_s_dot_w = dataAtPts->getFTensorSmallWL2Dot(nb_integration_pts);
786 auto w_l2_dot_dot_at_pts = dataAtPts->getSmallWL2DotDotAtPts();
787 const bool reset_w_l2_dot_dot =
788 w_l2_dot_dot_at_pts->size1() != nb_integration_pts ||
789 w_l2_dot_dot_at_pts->size2() != 3;
791 *w_l2_dot_dot_at_pts, nb_integration_pts);
792 if (reset_w_l2_dot_dot) {
793 w_l2_dot_dot_at_pts->clear();
795 auto t_s_dot_dot_w = dataAtPts->getFTensorSmallWL2DotDot(nb_integration_pts);
797 auto piola_scale = dataAtPts->piolaScale;
798 auto alpha_w = alphaW / piola_scale;
799 auto alpha_rho = alphaRho / piola_scale;
801 int nb_base_functions = data.
getN().size2();
805 auto get_ftensor1 = [](
auto &
v) {
817 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
819 auto t_nf = get_ftensor1(nF);
821 for (; bb != nb_dofs / 3; ++bb) {
822 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
823 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
824 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
828 for (; bb != nb_base_functions; ++bb)
839 int nb_integration_pts = getGaussPts().size2();
840 auto v = getVolume();
841 auto t_w = getFTensor0IntegrationWeight();
842 auto t_levi_kirchhoff =
843 dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts);
844 auto t_omega_grad_dot =
845 dataAtPts->getFTensorRotAxisGradDot(nb_integration_pts);
846 int nb_base_functions = data.
getN().size2();
852 auto get_ftensor1 = [](
auto &
v) {
858 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
861 auto t_nf = get_ftensor1(nF);
863 for (; bb != nb_dofs / 3; ++bb) {
864 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
865 t_nf(
k) += (
a * alphaOmega ) *
866 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
871 for (; bb != nb_base_functions; ++bb) {
885 int nb_integration_pts = data.
getN().size1();
886 auto v = getVolume();
887 auto t_w = getFTensor0IntegrationWeight();
889 int nb_base_functions = data.
getN().size2() / 3;
897 auto get_ftensor1 = [](
auto &
v) {
902 auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
904 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
906 auto t_nf = get_ftensor1(nF);
914 for (; bb != nb_dofs / 3; ++bb) {
915 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
920 for (; bb != nb_base_functions; ++bb)
933 int nb_integration_pts = data.
getN().size1();
934 auto v = getVolume();
935 auto t_w = getFTensor0IntegrationWeight();
937 int nb_base_functions = data.
getN().size2() / 9;
945 auto get_ftensor0 = [](
auto &
v) {
949 auto t_h = dataAtPts->getFTensorSmallH(nb_integration_pts);
951 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
953 auto t_nf = get_ftensor0(nF);
957 t_residuum(
i,
j) = t_h(
i,
j);
960 for (; bb != nb_dofs; ++bb) {
961 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
965 for (; bb != nb_base_functions; ++bb) {
978 int nb_integration_pts = data.
getN().size1();
979 auto v = getVolume();
980 auto t_w = getFTensor0IntegrationWeight();
981 auto t_w_l2 = dataAtPts->getFTensorSmallWL2(nb_integration_pts);
982 int nb_base_functions = data.
getN().size2() / 3;
985 auto get_ftensor1 = [](
auto &
v) {
990 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
992 auto t_nf = get_ftensor1(nF);
994 for (; bb != nb_dofs / 3; ++bb) {
995 double div_row_base = t_row_diff_base_fun(
i,
i);
996 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
998 ++t_row_diff_base_fun;
1000 for (; bb != nb_base_functions; ++bb) {
1001 ++t_row_diff_base_fun;
1015 int nb_integration_pts = getGaussPts().size2();
1018 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1020 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1021 if (tag_length != 9) {
1023 "Number of internal stress components should be 9 but is %d",
1028 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1029 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1030 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1031 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1033 auto get_internal_stress =
1035 dataAtPts->internalStressAtPts, nb_integration_pts);
1036 dataAtPts->internalStressAtPts.clear();
1037 auto t_internal_stress = get_internal_stress();
1040 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1041 t_internal_stress(
L) = t_const_stress(
L);
1042 ++t_internal_stress;
1053 int nb_integration_pts = getGaussPts().size2();
1056 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1058 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1059 if (tag_length != 9) {
1061 "Number of internal stress components should be 9 but is %d",
1065 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1068 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1071 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1074 auto get_internal_stress =
1076 dataAtPts->internalStressAtPts, nb_integration_pts);
1077 dataAtPts->internalStressAtPts.clear();
1078 auto t_internal_stress = get_internal_stress();
1083 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1084 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1085 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1086 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1090 ++t_internal_stress;
1102 int nb_integration_pts = data.
getN().size1();
1103 auto v = getVolume();
1104 auto t_w = getFTensor0IntegrationWeight();
1109 auto get_ftensor2 = [](
auto &
v) {
1111 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1114 auto t_internal_stress =
1115 dataAtPts->getFTensorInternalStress(nb_integration_pts);
1119 : getFEMethod()->ts_t;
1122 double scale = scalingMethodPtr->getScale(time);
1127 int nb_base_functions = data.
getN().size2();
1129 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1131 auto t_nf = get_ftensor2(nF);
1134 t_symm_stress(
i,
j) =
1135 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1138 t_residual(
L) = t_L(
i,
j,
L) * (
scale * t_symm_stress(
i,
j));
1141 for (; bb != nb_dofs / 6; ++bb) {
1142 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1146 for (; bb != nb_base_functions; ++bb)
1150 ++t_internal_stress;
1160 int nb_integration_pts = data.
getN().size1();
1161 auto v = getVolume();
1162 auto t_w = getFTensor0IntegrationWeight();
1164 auto get_ftensor2 = [](
auto &
v) {
1166 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1169 auto t_internal_stress =
1170 dataAtPts->getFTensorInternalStressVec(nb_integration_pts);
1179 : getFEMethod()->ts_t;
1182 double scale = scalingMethodPtr->getScale(time);
1184 int nb_base_functions = data.
getN().size2();
1186 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1188 auto t_nf = get_ftensor2(nF);
1191 t_residual(
L) = t_L(
M,
L) * (
scale * t_internal_stress(
M));
1194 for (; bb != nb_dofs / 6; ++bb) {
1195 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1199 for (; bb != nb_base_functions; ++bb)
1203 ++t_internal_stress;
1208template <AssemblyType A>
1214 for (
auto &bc : (*bcDispPtr)) {
1216 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1219 int nb_integration_pts = OP::getGaussPts().size2();
1220 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1221 auto t_w = OP::getFTensor0IntegrationWeight();
1222 int nb_base_functions = data.
getN().size2() / 3;
1229 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1231 scale *= scalingMethodsMap.at(bc.blockName)
1234 scale *= scalingMethodsMap.at(bc.blockName)
1235 ->getScale(OP::getFEMethod()->ts_t);
1239 <<
"No scaling method found for " << bc.blockName;
1246 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1247 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1249 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1251 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1255 for (; bb != nb_base_functions; ++bb)
1267 return OP::iNtegrate(data);
1275 for (
auto &bc : (*bcDispPtr)) {
1277 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1280 int nb_integration_pts = OP::getGaussPts().size2();
1281 auto t_w = OP::getFTensor0IntegrationWeight();
1282 int nb_base_functions = data.
getN().size2();
1285 if (!this->sourceVec) {
1287 "Source vector for OpTauStabilizationDispRhsBc is not set");
1289 if (data.
getN().size1() != nb_integration_pts) {
1291 "Number of integration points in data should be %d but is %d",
1292 nb_integration_pts, (
int)data.
getN().size1());
1294 if (nb_base_functions < nb_dofs /
SPACE_DIM) {
1296 "Number of base functions in data should be %d but is %d",
1304 *this->sourceVec, nb_integration_pts)();
1310 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1312 scale *= scalingMethodsMap.at(bc.blockName)
1315 scale *= scalingMethodsMap.at(bc.blockName)
1316 ->getScale(OP::getFEMethod()->ts_t);
1320 <<
"No scaling method found for " << bc.blockName;
1327 auto area = getMeasure();
1328 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1329 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1331 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1332 auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
1334 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1335 for (
auto ii = 0; ii !=
SPACE_DIM; ++ii) {
1337 t_nf(ii) += (tau_scale * t_row_base_fun) *
1338 (t_disp_val(ii) - t_bc_disp(ii));
1344 for (; bb != nb_base_functions; ++bb)
1363 for (
auto &bc : (*bcDispPtr)) {
1365 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1368 int nb_integration_pts = OP::getGaussPts().size2();
1369 auto t_w = OP::getFTensor0IntegrationWeight();
1370 int nb_base_functions = row_data.
getN().size2();
1376 auto get_t_vec = [&](
const int rr) {
1377 std::array<double *, SPACE_DIM> ptrs;
1379 ptrs[
i] = &OP::locMat(rr +
i,
i);
1384 auto area = getMeasure();
1385 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1386 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1388 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1390 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1393 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
1394 for (
int ii = 0; ii !=
SPACE_DIM; ++ii) {
1396 t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
1404 for (; rr != nb_base_functions; ++rr)
1420 for (
auto &bc : (*bcDispPtr)) {
1421 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1424 auto &v_analytical_expr = std::get<1>(analytical_data);
1427 int nb_integration_pts = OP::getGaussPts().size2();
1428 auto t_w = OP::getFTensor0IntegrationWeight();
1429 int nb_base_functions = data.
getN().size2();
1433 if (!this->sourceVec) {
1435 "Source vector for OpTauStabilizationOpAnalyticalDispBc is not "
1438 if (data.
getN().size1() != nb_integration_pts) {
1440 "Number of integration points in data should be %d but is %d",
1441 nb_integration_pts, (
int)data.
getN().size1());
1443 if (nb_base_functions < nb_dofs /
SPACE_DIM) {
1445 "Number of base functions in data should be at least %d but is "
1447 nb_dofs /
SPACE_DIM, nb_base_functions);
1453 *this->sourceVec, nb_integration_pts)();
1456 auto area = getMeasure();
1457 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1458 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1460 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1461 auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
1463 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1464 for (
auto ii = 0; ii !=
SPACE_DIM; ++ii) {
1467 (tau_scale * t_row_base_fun) * (t_disp_val(ii) - t_bc_disp(ii));
1473 for (; bb != nb_base_functions; ++bb)
1492 for (
auto &bc : (*bcDispPtr)) {
1493 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1495 int nb_integration_pts = OP::getGaussPts().size2();
1496 auto t_w = OP::getFTensor0IntegrationWeight();
1497 int nb_base_functions = row_data.
getN().size2();
1500 auto get_t_vec = [&](
const int rr) {
1501 std::array<double *, SPACE_DIM> ptrs;
1503 ptrs[
i] = &OP::locMat(rr +
i,
i);
1508 auto area = getMeasure();
1509 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1510 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1512 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1514 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1517 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
1518 for (
int ii = 0; ii !=
SPACE_DIM; ++ii) {
1520 t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
1528 for (; rr != nb_base_functions; ++rr)
1540template <AssemblyType A>
1548 double time = OP::getFEMethod()->ts_t;
1556 for (
auto &bc : (*bcRotPtr)) {
1558 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1560 int nb_integration_pts = OP::getGaussPts().size2();
1561 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1562 auto t_w = OP::getFTensor0IntegrationWeight();
1564 int nb_base_functions = data.
getN().size2() / 3;
1567 auto get_ftensor1 = [](
auto &
v) {
1580 auto get_rotation_angle = [&]() {
1581 double theta = bc.theta;
1582 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1583 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1588 auto get_rotation = [&](
auto theta) {
1590 if (bc.vals.size() == 7) {
1591 t_omega(0) = bc.vals[4];
1592 t_omega(1) = bc.vals[5];
1593 t_omega(2) = bc.vals[6];
1596 t_omega(
i) = OP::getFTensor1Normal()(
i);
1598 if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
1602 <<
"Rotation axis is zero vector for block " << bc.blockName
1603 <<
". This may lead to unexpected results.";
1605 t_omega(
i) *= theta;
1612 auto t_R = get_rotation(get_rotation_angle());
1613 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1615 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1617 t_delta(
i) = t_center(
i) - t_coords(
i);
1619 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1621 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1623 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1624 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1628 for (; bb != nb_base_functions; ++bb)
1641 return OP::iNtegrate(data);
1650 double time = OP::getFEMethod()->ts_t;
1658 for (
auto &bc : (*bcRotPtr)) {
1660 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1662 int nb_integration_pts = OP::getGaussPts().size2();
1663 auto t_w = OP::getFTensor0IntegrationWeight();
1665 int nb_base_functions = data.
getN().size2();
1668 auto get_ftensor1 = [](
auto &
v) {
1681 auto get_rotation_angle = [&]() {
1682 double theta = bc.theta;
1683 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1684 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1689 auto get_rotation = [&](
auto theta) {
1691 if (bc.vals.size() == 7) {
1692 t_omega(0) = bc.vals[4];
1693 t_omega(1) = bc.vals[5];
1694 t_omega(2) = bc.vals[6];
1697 t_omega(
i) = OP::getFTensor1Normal()(
i);
1699 if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
1702 t_omega(
i) *= theta;
1709 auto area = getMeasure();
1710 auto t_R = get_rotation(get_rotation_angle());
1711 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1714 *this->sourceVec, nb_integration_pts)();
1716 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1718 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1721 t_delta(
i) = t_center(
i) - t_coords(
i);
1723 t_bc_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1725 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1727 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1729 (tau_scale * t_row_base_fun) * (t_disp_val(
i) - t_bc_disp(
i));
1733 for (; bb != nb_base_functions; ++bb)
1752 for (
auto &bc : (*bcRotPtr)) {
1754 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1757 int nb_integration_pts = OP::getGaussPts().size2();
1758 auto t_w = OP::getFTensor0IntegrationWeight();
1759 int nb_base_functions = row_data.
getN().size2();
1765 auto get_t_vec = [&](
const int rr) {
1766 std::array<double *, SPACE_DIM> ptrs;
1768 ptrs[
i] = &OP::locMat(rr +
i,
i);
1773 auto area = getMeasure();
1774 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1775 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1777 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1779 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1782 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
1783 for (
int ii = 0; ii !=
SPACE_DIM; ++ii) {
1784 t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
1791 for (; rr != nb_base_functions; ++rr)
1803template <AssemblyType A>
1807 double time = OP::getFEMethod()->ts_t;
1815 for (
auto &bc : (*bcDispPtr)) {
1817 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1819 for (
auto &bd : (*brokenBaseSideDataPtr)) {
1823 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1824 auto t_w = OP::getFTensor0IntegrationWeight();
1832 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1833 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1836 <<
"No scaling method found for " << bc.blockName;
1840 double val =
scale * bc.val;
1843 int nb_integration_pts = OP::getGaussPts().size2();
1844 int nb_base_functions = data.
getN().size2();
1846 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1849 t_N(
i) = t_normal(
i);
1853 t_P(
i,
j) = t_N(
i) * t_N(
j);
1858 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
1862 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
1864 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1866 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1867 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
1871 for (; bb != nb_base_functions; ++bb)
1885template <AssemblyType A>
1891 double time = OP::getFEMethod()->ts_t;
1896 int row_nb_dofs = row_data.
getIndices().size();
1897 int col_nb_dofs = col_data.
getIndices().size();
1898 auto &locMat = OP::locMat;
1899 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1905 for (
auto &bc : (*bcDispPtr)) {
1907 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1909 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1910 auto t_w = OP::getFTensor0IntegrationWeight();
1916 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1917 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1920 <<
"No scaling method found for " << bc.blockName;
1923 int nb_integration_pts = OP::getGaussPts().size2();
1924 int row_nb_dofs = row_data.
getIndices().size();
1925 int col_nb_dofs = col_data.
getIndices().size();
1926 int nb_base_functions = row_data.
getN().size2();
1931 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1934 t_N(
i) = t_normal(
i);
1938 t_P(
i,
j) = t_N(
i) * t_N(
j);
1941 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
1944 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
1945 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1948 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
1949 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
1956 for (; rr != nb_base_functions; ++rr)
1963 locMat *= OP::getMeasure();
1969template <AssemblyType A>
1975 double time = OP::getFEMethod()->ts_t;
1980 int row_nb_dofs = row_data.
getIndices().size();
1981 int col_nb_dofs = col_data.
getIndices().size();
1982 auto &locMat = OP::locMat;
1983 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
1989 for (
auto &bc : (*bcDispPtr)) {
1991 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1993 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1994 auto t_w = OP::getFTensor0IntegrationWeight();
2003 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2004 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2007 <<
"No scaling method found for " << bc.blockName;
2010 int nb_integration_pts = OP::getGaussPts().size2();
2011 int nb_base_functions = row_data.
getN().size2();
2014 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2017 t_N(
i) = t_normal(
i);
2021 t_P(
i,
j) = t_N(
i) * t_N(
j);
2026 t_d_res(
i,
j) = t_Q(
i,
j);
2029 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2030 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2033 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2035 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
2042 for (; rr != nb_base_functions; ++rr)
2049 locMat *= OP::getMeasure();
2056 return OP::iNtegrate(data);
2061 return OP::iNtegrate(row_data, col_data);
2066 return OP::iNtegrate(row_data, col_data);
2069template <AssemblyType A>
2073 double time = OP::getFEMethod()->ts_t;
2081 for (
auto &bc : (*bcDispPtr)) {
2083 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2088 auto [block_name, v_analytical_expr] =
2093 int nb_integration_pts = OP::getGaussPts().size2();
2094 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2095 auto t_w = OP::getFTensor0IntegrationWeight();
2096 int nb_base_functions = data.
getN().size2() / 3;
2098 auto t_coord = OP::getFTensor1CoordsAtGaussPts();
2106 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2107 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
2110 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
2112 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
2116 for (; bb != nb_base_functions; ++bb)
2130 return OP::iNtegrate(data);
2139 int nb_integration_pts = getGaussPts().size2();
2140 int nb_base_functions = data.
getN().size2();
2142 double time = getFEMethod()->ts_t;
2148 if (this->locF.size() != nb_dofs)
2150 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2153 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2156 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
2158 auto t_w = getFTensor0IntegrationWeight();
2159 auto t_coords = getFTensor1CoordsAtGaussPts();
2160 auto t_normal = getFTensor1NormalsAtGaussPts();
2162 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2164 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2166 double a = sqrt(t_normal(
i) * t_normal(
i));
2168 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2169 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2171 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2173 (time_scale *
a * t_w * t_row_base * tau) * (t_val(
i) *
scale);
2178 for (; rr != nb_base_functions; ++rr)
2189 for (
auto &bc : *(bcData)) {
2190 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2192 double time_scale = 1;
2193 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2194 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2200 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
2204 return -y * (y - 1) / 0.25;
2206 CHKERR integrate_rhs(bc, calc_tau, time_scale);
2209 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2223 int nb_integration_pts = getGaussPts().size2();
2224 int nb_base_functions = data.
getN().size2();
2226 double time = getFEMethod()->ts_t;
2232 if (this->locF.size() != nb_dofs)
2234 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2237 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2242 auto t_w = getFTensor0IntegrationWeight();
2243 auto t_coords = getFTensor1CoordsAtGaussPts();
2244 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2245 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2249 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2251 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2263 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2266 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2267 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2268 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2270 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2272 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2274 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2276 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2277 t_f(
i) += t_row_base * t_val(
i);
2282 for (; rr != nb_base_functions; ++rr)
2297 for (
auto &bc : *(bcData)) {
2298 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2300 double time_scale = 1;
2301 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2302 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2308 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2315template <AssemblyType A>
2326 double time = OP::getFEMethod()->ts_t;
2331 int nb_base_functions = row_data.
getN().size2();
2332 int row_nb_dofs = row_data.
getIndices().size();
2333 int col_nb_dofs = col_data.
getIndices().size();
2334 int nb_integration_pts = OP::getGaussPts().size2();
2335 auto &locMat = OP::locMat;
2336 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2339 auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2344 auto t_w = OP::getFTensor0IntegrationWeight();
2345 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2346 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2347 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2352 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2362 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2363 auto t_val = time_scale * t_w * tau * val;
2366 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2367 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2370 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2372 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2373 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2374 t_kd(
j,
l) * t_diff_col_base(N0)
2378 (FTensor::levi_civita<double>(
i,
j,
k) *
2379 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2380 t_kd(
k,
l) * t_diff_col_base(N1);
2382 t_mat(
i,
j) += (t_w * t_row_base) * t_val * t_normal_du(
i,
j);
2389 for (; rr != nb_base_functions; ++rr)
2405 for (
auto &bc : *(bcData)) {
2406 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2408 double time_scale = 1;
2409 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2410 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2416 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2426 return OP::iNtegrate(row_data, col_data);
2435 int nb_integration_pts = getGaussPts().size2();
2436 int nb_base_functions = data.
getN().size2();
2438 double time = getFEMethod()->ts_t;
2444 if (this->locF.size() != nb_dofs)
2446 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2451 for (
auto &bc : *(bcData)) {
2452 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2456 auto [block_name, v_analytical_expr] =
2460 auto t_w = getFTensor0IntegrationWeight();
2461 auto t_coords = getFTensor1CoordsAtGaussPts();
2463 double scale = (piolaScalePtr) ? 1. / (*piolaScalePtr) : 1.0;
2465 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2467 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2469 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2470 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2475 for (; rr != nb_base_functions; ++rr)
2481 this->locF *= getMeasure();
2490 int nb_integration_pts = row_data.
getN().size1();
2491 int row_nb_dofs = row_data.
getIndices().size();
2492 int col_nb_dofs = col_data.
getIndices().size();
2493 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2495 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2498 auto v = getVolume();
2499 auto t_w = getFTensor0IntegrationWeight();
2500 int row_nb_base_functions = row_data.
getN().size2();
2502 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2505 for (; rr != row_nb_dofs / 3; ++rr) {
2507 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2508 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2509 double div_col_base = t_col_diff_base_fun(
i,
i);
2510 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2512 ++t_col_diff_base_fun;
2516 for (; rr != row_nb_base_functions; ++rr)
2527 if (alphaW < std::numeric_limits<double>::epsilon() &&
2528 alphaRho < std::numeric_limits<double>::epsilon())
2531 const int nb_integration_pts = row_data.
getN().size1();
2532 const int row_nb_dofs = row_data.
getIndices().size();
2533 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2535 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
2541 auto v = getVolume();
2542 auto t_w = getFTensor0IntegrationWeight();
2544 auto piola_scale = dataAtPts->piolaScale;
2545 auto alpha_w = alphaW / piola_scale;
2546 auto alpha_rho = alphaRho / piola_scale;
2548 int row_nb_base_functions = row_data.
getN().size2();
2551 double ts_scale = alpha_w * getTSa();
2552 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
2553 ts_scale += alpha_rho * getTSaa();
2555 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2556 double a =
v * t_w * ts_scale;
2559 for (; rr != row_nb_dofs / 3; ++rr) {
2562 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2563 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
2564 const double b =
a * t_row_base_fun * t_col_base_fun;
2573 for (; rr != row_nb_base_functions; ++rr)
2592 int nb_integration_pts = row_data.
getN().size1();
2593 int row_nb_dofs = row_data.
getIndices().size();
2594 int col_nb_dofs = col_data.
getIndices().size();
2595 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2598 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2600 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2602 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2604 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2606 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2608 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
2611 auto v = getVolume();
2612 auto t_w = getFTensor0IntegrationWeight();
2614 int row_nb_base_functions = row_data.
getN().size2();
2617 auto t_approx_P_adjoint_log_du_dP =
2618 dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
2620 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2623 for (; rr != row_nb_dofs / 6; ++rr) {
2626 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2628 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2630 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
2638 for (; rr != row_nb_base_functions; ++rr)
2641 ++t_approx_P_adjoint_log_du_dP;
2657 int nb_integration_pts = row_data.
getN().size1();
2658 int row_nb_dofs = row_data.
getIndices().size();
2659 int col_nb_dofs = col_data.
getIndices().size();
2660 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2662 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
2666 auto v = getVolume();
2667 auto t_w = getFTensor0IntegrationWeight();
2670 int row_nb_base_functions = row_data.
getN().size2();
2672 auto t_approx_P_adjoint_log_du_dP =
2673 dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
2675 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2678 for (; rr != row_nb_dofs / 6; ++rr) {
2679 auto t_m = get_ftensor2(
K, 6 * rr, 0);
2680 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2681 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2683 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
2690 for (; rr != row_nb_base_functions; ++rr)
2693 ++t_approx_P_adjoint_log_du_dP;
2705 int nb_integration_pts = getGaussPts().size2();
2706 int row_nb_dofs = row_data.
getIndices().size();
2707 int col_nb_dofs = col_data.
getIndices().size();
2708 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
2711 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2713 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2715 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
2717 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
2719 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
2721 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
2731 auto v = getVolume();
2732 auto t_w = getFTensor0IntegrationWeight();
2733 auto t_approx_P_adjoint_log_du_domega =
2734 dataAtPts->getFTensorAdjointPdUdOmega(nb_integration_pts);
2736 int row_nb_base_functions = row_data.
getN().size2();
2739 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2743 for (; rr != row_nb_dofs / 6; ++rr) {
2745 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2746 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2747 double v =
a * t_row_base_fun * t_col_base_fun;
2748 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
2755 for (; rr != row_nb_base_functions; ++rr)
2759 ++t_approx_P_adjoint_log_du_domega;
2768 int nb_integration_pts = getGaussPts().size2();
2769 int row_nb_dofs = row_data.
getIndices().size();
2770 int col_nb_dofs = col_data.
getIndices().size();
2771 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2775 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
2776 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
2778 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
2779 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
2781 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
2782 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
2790 auto v = getVolume();
2791 auto t_w = getFTensor0IntegrationWeight();
2792 auto t_levi_kirchhoff_du =
2793 dataAtPts->getFTensorLeviKirchhoffdLogStretch(nb_integration_pts);
2794 int row_nb_base_functions = row_data.
getN().size2();
2796 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2799 for (; rr != row_nb_dofs / 3; ++rr) {
2800 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2801 const double b =
a * t_row_base_fun;
2803 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
2804 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
2810 for (; rr != row_nb_base_functions; ++rr) {
2814 ++t_levi_kirchhoff_du;
2830 int nb_integration_pts = getGaussPts().size2();
2831 int row_nb_dofs = row_data.
getIndices().size();
2832 int col_nb_dofs = col_data.
getIndices().size();
2833 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2836 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2838 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2840 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2845 auto v = getVolume();
2846 auto t_w = getFTensor0IntegrationWeight();
2848 int row_nb_base_functions = row_data.
getN().size2();
2850 auto t_levi_kirchhoff_dP =
2851 dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
2853 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2856 for (; rr != row_nb_dofs / 3; ++rr) {
2857 double b =
a * t_row_base_fun;
2859 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2860 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2861 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
2867 for (; rr != row_nb_base_functions; ++rr) {
2872 ++t_levi_kirchhoff_dP;
2880 int nb_integration_pts = getGaussPts().size2();
2881 int row_nb_dofs = row_data.
getIndices().size();
2882 int col_nb_dofs = col_data.
getIndices().size();
2884 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2886 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2893 auto v = getVolume();
2894 auto t_w = getFTensor0IntegrationWeight();
2895 auto t_levi_kirchoff_dP =
2896 dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
2898 int row_nb_base_functions = row_data.
getN().size2();
2901 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2904 for (; rr != row_nb_dofs / 3; ++rr) {
2905 double b =
a * t_row_base_fun;
2906 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2907 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2908 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2909 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
2916 for (; rr != row_nb_base_functions; ++rr) {
2920 ++t_levi_kirchoff_dP;
2928 int nb_integration_pts = getGaussPts().size2();
2929 int row_nb_dofs = row_data.
getIndices().size();
2930 int col_nb_dofs = col_data.
getIndices().size();
2931 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
2934 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2936 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2938 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2951 auto v = getVolume();
2952 auto ts_a = std::abs(alphaOmega) > std::numeric_limits<double>::epsilon()
2955 auto t_w = getFTensor0IntegrationWeight();
2956 auto t_levi_kirchhoff_domega =
2957 dataAtPts->getFTensorLeviKirchhoffdOmega(nb_integration_pts);
2958 int row_nb_base_functions = row_data.
getN().size2();
2963 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2965 double c = (
a * alphaOmega) * (ts_a );
2968 for (; rr != row_nb_dofs / 3; ++rr) {
2969 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2970 const double b =
a * t_row_base_fun;
2973 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2974 t_m(
k,
l) -= (b * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
2975 t_m(
k,
l) +=
t_kd(
k,
l) * (
c * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
2983 for (; rr != row_nb_base_functions; ++rr) {
2988 ++t_levi_kirchhoff_domega;
2996 if (hasNonhomogeneousMatBlock) {
2998 integrateImpl<size_symm * size_symm>(row_data, col_data));
3010 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3013 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3015 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3017 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3022 int nb_integration_pts = getGaussPts().size2();
3023 int row_nb_dofs = row_data.
getIndices().size();
3024 int col_nb_dofs = col_data.
getIndices().size();
3026 auto v = getVolume();
3027 auto t_w = getFTensor0IntegrationWeight();
3028 int row_nb_base_functions = row_data.
getN().size2() / 3;
3036 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
3039 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3043 for (; rr != row_nb_dofs / 3; ++rr) {
3045 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3046 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3047 t_m(
i,
k) -=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
3055 for (; rr != row_nb_base_functions; ++rr)
3068 if (hasNonhomogeneousMatBlock) {
3070 integrateImpl<size_symm * size_symm>(row_data, col_data));
3083 int nb_integration_pts = getGaussPts().size2();
3084 int row_nb_dofs = row_data.
getIndices().size();
3085 int col_nb_dofs = col_data.
getIndices().size();
3087 auto v = getVolume();
3088 auto t_w = getFTensor0IntegrationWeight();
3089 int row_nb_base_functions = row_data.
getN().size2() / 9;
3097 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
3100 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3104 for (; rr != row_nb_dofs; ++rr) {
3106 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
3108 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
3115 for (; rr != row_nb_base_functions; ++rr)
3126 if (hasNonhomogeneousMatBlock) {
3128 integrateImpl<size_symm * size_symm>(row_data, col_data));
3141 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3144 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
3149 int nb_integration_pts = getGaussPts().size2();
3150 int row_nb_dofs = row_data.
getIndices().size();
3151 int col_nb_dofs = col_data.
getIndices().size();
3153 auto v = getVolume();
3154 auto t_w = getFTensor0IntegrationWeight();
3155 int row_nb_base_functions = row_data.
getN().size2() / 9;
3165 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(dataAtPts->matInvD);
3168 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3171 auto t_m = get_ftensor1(
K, 0, 0);
3174 for (; rr != row_nb_dofs; ++rr) {
3176 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3177 t_m(
k) -=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
3185 for (; rr != row_nb_base_functions; ++rr)
3197 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3200 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
3205 int nb_integration_pts = getGaussPts().size2();
3206 int row_nb_dofs = row_data.
getIndices().size();
3207 int col_nb_dofs = col_data.
getIndices().size();
3209 auto v = getVolume();
3210 auto t_w = getFTensor0IntegrationWeight();
3211 int row_nb_base_functions = row_data.
getN().size2() / 9;
3217 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3220 auto t_m = get_ftensor1(
K, 0, 0);
3223 for (; rr != row_nb_dofs; ++rr) {
3225 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3226 t_m(
k) +=
a * t_row_base(
k,
l) * t_col_base(
l);
3234 for (; rr != row_nb_base_functions; ++rr)
3253 int nb_integration_pts = row_data.
getN().size1();
3254 int row_nb_dofs = row_data.
getIndices().size();
3255 int col_nb_dofs = col_data.
getIndices().size();
3257 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3260 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3262 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3264 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3269 auto v = getVolume();
3270 auto t_w = getFTensor0IntegrationWeight();
3271 int row_nb_base_functions = row_data.
getN().size2() / 3;
3274 auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
3276 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3280 for (; rr != row_nb_dofs / 3; ++rr) {
3283 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
3286 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3287 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3288 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
3296 for (; rr != row_nb_base_functions; ++rr)
3316 int nb_integration_pts = row_data.
getN().size1();
3317 int row_nb_dofs = row_data.
getIndices().size();
3318 int col_nb_dofs = col_data.
getIndices().size();
3320 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3322 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3325 auto v = getVolume();
3326 auto t_w = getFTensor0IntegrationWeight();
3327 int row_nb_base_functions = row_data.
getN().size2() / 9;
3330 auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
3331 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3335 for (; rr != row_nb_dofs; ++rr) {
3338 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3341 auto t_m = get_ftensor2(
K, rr, 0);
3342 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3343 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3351 for (; rr != row_nb_base_functions; ++rr)
3364 if (tagSense != getSkeletonSense())
3367 auto get_tag = [&](
auto name) {
3368 auto &mob = getPtrFE()->mField.get_moab();
3374 auto get_tag_value = [&](
auto &&tag,
int dim) {
3375 auto &mob = getPtrFE()->mField.get_moab();
3376 auto face = getSidePtrFE()->getFEEntityHandle();
3377 std::vector<double> value(dim);
3378 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3382 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3383 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3385 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
3386 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3391 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3392 Tag th_detF = create_tag(
"detF", 1);
3393 Tag th_traction = create_tag(
"traction", 3);
3394 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3396 Tag th_energy = create_tag(
"Energy", 1);
3397 Tag th_young_modulus = create_tag(
"YoungModulus", 1);
3399 const auto nb_gauss_pts = getGaussPts().size2();
3400 auto t_w = dataAtPts->getFTensorSmallWL2(nb_gauss_pts);
3401 auto t_h = dataAtPts->getFTensorSmallH(nb_gauss_pts);
3402 auto t_approx_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
3404 auto t_normal = getFTensor1NormalsAtGaussPts();
3405 auto t_disp = dataAtPts->getFTensorSmallWH1(nb_gauss_pts);
3409 if (dataAtPts->energyAtPts.size() == 0) {
3411 dataAtPts->energyAtPts.resize(nb_gauss_pts);
3412 dataAtPts->energyAtPts.clear();
3432 auto set_float_precision = [](
const double x) {
3433 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3440 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3442 v = set_float_precision(
v);
3443 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
3450 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3453 for (
auto &
a :
v.data())
3454 a = set_float_precision(
a);
3455 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3456 &*
v.data().begin());
3464 &
m(0, 0), &
m(0, 1), &
m(0, 2),
3466 &
m(1, 0), &
m(1, 1), &
m(1, 2),
3468 &
m(2, 0), &
m(2, 1), &
m(2, 2));
3470 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
3472 t_m(
i,
j) = t_d(
i,
j);
3473 for (
auto &
v :
m.data())
3474 v = set_float_precision(
v);
3475 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
3476 &*
m.data().begin());
3480 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
3483 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
3485 t_traction(
i) *= tagSense;
3486 CHKERR save_vec_tag(th_traction, t_traction, gg);
3488 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
3489 if (!std::isfinite(u_error))
3491 CHKERR save_scal_tag(th_disp_error, u_error, gg);
3492 CHKERR save_scal_tag(th_energy, t_energy, gg);
3494 CHKERR save_scal_tag(th_young_modulus, t_youngs_modulus, gg);
3498 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
3499 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
3500 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
3509 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3510 std::vector<FieldSpace> spaces, std::string geom_field_name,
3511 boost::shared_ptr<Range> crack_front_edges_ptr) {
3514 constexpr bool scale_l2 =
false;
3518 "Scale L2 Ainsworth Legendre base is not implemented");
3527 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3528 std::vector<FieldSpace> spaces, std::string geom_field_name,
3529 boost::shared_ptr<Range> crack_front_edges_ptr) {
3532 constexpr bool scale_l2 =
false;
3536 "Scale L2 Ainsworth Legendre base is not implemented");
3545 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
3546 std::vector<FieldSpace> spaces, std::string geom_field_name,
3547 boost::shared_ptr<Range> crack_front_edges_ptr,
3548 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
3549 boost::shared_ptr<MatrixDouble> inv_jac) {
3552 if (!geom_field_name.empty()) {
3554 auto jac = boost::make_shared<MatrixDouble>();
3555 auto det = boost::make_shared<VectorDouble>();
3558 geom_field_name, jac));
3565 constexpr bool scale_l2_ainsworth_legendre_base =
false;
3567 if (scale_l2_ainsworth_legendre_base) {
3575 boost::shared_ptr<MatrixDouble> jac,
3576 boost::shared_ptr<Range> edges_ptr)
3585 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
3588 return OP::doWork(side,
type, data);
3593 boost::shared_ptr<Range> edgesPtr;
3596 if (!geom_field_name.empty()) {
3597 auto jac = boost::make_shared<MatrixDouble>();
3598 auto det = boost::make_shared<VectorDouble>();
3600 geom_field_name, jac,
3602 : boost::make_shared<
Range>()));
3677 const auto nb_gauss_pts = getGaussPts().size2();
3679 dataAtPts->faceMaterialForceAtPts, nb_gauss_pts);
3680 dataAtPts->normalPressureAtPts.resize(nb_gauss_pts,
false);
3681 if (getNinTheLoop() == 0) {
3682 dataAtPts->faceMaterialForceAtPts.clear();
3683 dataAtPts->normalPressureAtPts.clear();
3685 auto loop_size = getLoopSize();
3686 if (loop_size == 1) {
3687 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
3688 auto pstatus = numebered_fe_ptr->getPStatus();
3689 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
3696 auto t_normal = getFTensor1NormalsAtGaussPts();
3697 auto t_T = dataAtPts->getFTensorFaceMaterialForce(
3701 auto t_P = dataAtPts->getFTensorApproxP(nb_gauss_pts);
3702 auto t_u_gamma = dataAtPts->getFTensorSmallHybridDisp(nb_gauss_pts);
3703 auto t_grad_u_gamma = dataAtPts->getFTensorGradHybridDisp(nb_gauss_pts);
3704 auto t_strain = dataAtPts->getFTensorLogStretch(nb_gauss_pts);
3705 auto t_omega = dataAtPts->getFTensorRotAxis(nb_gauss_pts);
3727 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3728 t_N(
I) = t_normal(
I);
3731 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
3733 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
3735 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3738 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3741 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3748 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3751 t_N(
I) = t_normal(
I);
3756 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
3760 t_grad_u_gamma(
i,
J) +
3761 (2 * t_R(
i,
K) * t_N(
K) - (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
3764 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
3767 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
3771 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
3780 "Grffith energy release "
3781 "selector not implemented");
3785 auto side_fe_ptr = getSidePtrFE();
3786 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
3787 auto pstatus = side_fe_mi_ptr->getPStatus();
3789 auto owner = side_fe_mi_ptr->getOwnerProc();
3791 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
3792 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
3793 << getNinTheLoop() <<
" loop size " << getLoopSize();
3805 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
3806 auto pstatus = fe_mi_ptr->getPStatus();
3808 auto owner = fe_mi_ptr->getOwnerProc();
3810 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
3811 <<
" " << getPtrFE()->mField.get_comm_rank();
3819 double face_pressure = 0.;
3820 auto t_T = dataAtPts->getFTensorFaceMaterialForce(
3821 getGaussPts().size2());
3824 auto t_w = getFTensor0IntegrationWeight();
3825 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
3826 t_face_T(
I) += t_w * t_T(
I);
3827 face_pressure += t_w * t_p;
3832 t_face_T(
I) *= getMeasure();
3833 face_pressure *= getMeasure();
3835 auto get_tag = [&](
auto name,
auto dim) {
3836 auto &moab = getPtrFE()->mField.get_moab();
3838 double def_val[] = {0., 0., 0.};
3839 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
3840 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
3845 auto set_tag = [&](
auto &&tag,
auto ptr) {
3846 auto &moab = getPtrFE()->mField.get_moab();
3847 auto face = getPtrFE()->getFEEntityHandle();
3848 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
3851 set_tag(
get_tag(
"MaterialForce", 3), &t_face_T(0));
3852 set_tag(
get_tag(
"FacePressure", 1), &face_pressure);
3857template <
typename OP_PTR>
3858std::tuple<std::string, MatrixDouble>
3860 const std::string block_name) {
3862 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
3864 auto ts_time = op_ptr->getTStime();
3865 auto ts_time_step = op_ptr->getTStimeStep();
3872 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
3873 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
3875 auto v_analytical_expr =
3877 m_ref_coords, m_ref_normals, block_name);
3879 if (PetscUnlikely(!v_analytical_expr.size2())) {
3881 "Analytical expression is empty or does not exist, "
3882 "check python file");
3885 return std::make_tuple(block_name, v_analytical_expr);
3889 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3890 boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
3891 boost::shared_ptr<Range> ents_ptr)
3892 :
OP(broken_base_side_data, ents_ptr) {
3893 this->sourceVec = vec;
3894 this->betaCoeff = beta_coeff;
3898 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3899 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3900 :
OP(broken_base_side_data, ents_ptr) {
3901 this->sourceVec = boost::shared_ptr<MatrixDouble>();
3902 this->betaCoeff = beta_coeff;
3911 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
3916 if (!brokenBaseSideData) {
3921 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
3929 OP::nbIntegrationPts = OP::getGaussPts().size2();
3931 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
3933 OP::locF.resize(OP::nbRows,
false);
3936 CHKERR this->iNtegrate(row_data);
3938 CHKERR this->aSsemble(row_data);
3942 switch (OP::opType) {
3944 for (
auto &bd : *brokenBaseSideData) {
3946 boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
3947 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
3948 this->sourceVec.reset();
3953 (std::string(
"wrong op type ") +
3954 OpBaseDerivativesBase::OpTypeNames[OP::opType])
3962 const std::string row_field,
3963 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3964 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
3965 :
OP(row_field, boost::shared_ptr<
MatrixDouble>(), beta_coeff, ents_ptr),
3966 brokenBaseSideDataPtr(broken_base_side_data) {
3967 this->betaCoeff = beta_coeff;
3973 for (
auto &bd : (*brokenBaseSideDataPtr)) {
3978 if (this->sourceVec->size2() !=
SPACE_DIM) {
3980 "Inconsistent size of the source vector");
3982 if (this->sourceVec->size1() != OP::getGaussPts().size2()) {
3984 "Inconsistent size of the source vector");
3988 CHKERR OP::iNtegrate(data);
3990 this->sourceVec.reset();
3996 std::string row_field,
3997 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
3998 ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
3999 boost::shared_ptr<Range> ents_ptr)
4000 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
4002 this->betaCoeff = beta;
4007 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4008 ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
4009 :
OP(broken_base_side_data, ents_ptr) {
4011 this->betaCoeff = beta;
4012 OP::assembleTranspose =
false;
4013 OP::onlyTranspose =
false;
4022 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
4027 if (!brokenBaseSideData) {
4032 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
4033 EntityType col_type,
4038 auto check_if_assemble_transpose = [&] {
4040 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
4044 }
else if (OP::assembleTranspose) {
4050 OP::rowSide = row_side;
4051 OP::rowType = row_type;
4052 OP::colSide = col_side;
4053 OP::colType = col_type;
4055 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
4057 CHKERR this->iNtegrate(row_data, col_data);
4058 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
4062 switch (OP::opType) {
4065 for (
auto &bd : *brokenBaseSideData) {
4068 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
4070 "base functions not set");
4074 OP::nbRows = bd.getData().getIndices().size();
4077 OP::nbIntegrationPts = OP::getGaussPts().size2();
4078 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
4086 bd.getSide(), bd.getSide(),
4089 bd.getType(), bd.getType(),
4092 bd.getData(), bd.getData()
4101 (std::string(
"wrong op type ") +
4102 OpBaseDerivativesBase::OpTypeNames[OP::opType])
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
Tensor1< T, Tensor_Dim > normalize()
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_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)
static Tag get_tag(moab::Interface &moab, std::string tag_name, int size)
std::tuple< std::string, MatrixDouble > getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
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
implementation of Data Operators for Forces and Sources
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
decltype(GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
constexpr IntegrationType I
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static PetscBool l2UserBaseScale
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool physicalTimeFlg
static double currentPhysicalTime
static boost::function< double(const double)> f
static PetscBool setSingularity
static bool hasNonHomogeneousMaterialBlock
static boost::function< double(const double)> d_f
static enum EnergyReleaseSelector energyReleaseSelector
static boost::function< double(const double)> inv_f
static auto diffDiffExp(A &&t_w_vee, B &&theta)
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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
auto getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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.
auto getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Operator for inverting matrices at integration points.
Scale base functions by inverses of measure of element.
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseBrokenBase(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpBrokenBaseTimesBrokenDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
OpBrokenBaseTimesHybridDisp(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< MatrixDouble > vec, ScalarFun beta_coeff=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Caluclate face material force and normal pressure at gauss points.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpHybridBaseTimesBrokenDisp(const std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta_coeff, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< std::vector< BrokenBaseSideData > > brokenBaseSideDataPtr
OpHyrbridBaseBrokenBase(std::string row_field, boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, ScalarFun beta, const bool assmb_transpose, const bool only_transpose, boost::shared_ptr< Range > ents_ptr=nullptr)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
