14#include <boost/math/constants/constants.hpp>
27 boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
28 std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
30 externalPressurePtr(
std::move(external_pressure_ptr)),
31 externalStrainVecPtr(
std::move(external_strain_vec_ptr)),
32 scalingMethodsMap(
std::move(smv)) {
43 "External-pressure integration-point vector is null");
46 const int nb_integration_pts =
getGaussPts().size2();
63 const std::regex analytical_pattern(
"(.*)ANALYTICAL_EXTERNALSTRAIN(.*)");
66 if (block.ents.find(fe_ent) == block.ents.end()) {
69 if (!std::isfinite(block.bulkModulusK)) {
71 "External-strain bulk modulus is not finite in block %s",
72 block.blockName.c_str());
76 if (std::regex_match(block.blockName, analytical_pattern)) {
78 external_strain = analytical_externalstrain_function(
79 time_step, time, nb_integration_pts, reference_coordinates,
87 "Scaling method is null for external-strain block %s",
88 block.blockName.c_str());
90 scale = it->second->getScale(time);
93 <<
"No scaling method found for " << block.blockName;
95 std::fill(external_strain.begin(), external_strain.end(),
99 if (external_strain.size() != nb_integration_pts) {
101 "Wrong number of analytical external-strain integration points");
103 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
104 const double q = 3 * block.bulkModulusK * external_strain[gg];
105 if (!std::isfinite(q)) {
107 "External pressure q is not finite in block %s",
108 block.blockName.c_str());
110 (*externalPressurePtr)[gg] += q;
113 "Accumulated external pressure q is not finite in block %s",
114 block.blockName.c_str());
135 auto get_eshelby_stress =
137 DL>::size(
dataAtPts->SigmaAtPts, nb_integration_pts);
138 auto t_eshelby_stress = get_eshelby_stress();
142 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
143 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
153template <
typename TInvD,
typename TRotation,
typename TDiffRotation,
156 TDiffRotation &t_diff_R, TStress &t_P) {
159 t_d_rotated_p_d_omega;
160 t_d_rotated_p_d_omega(
i,
j,
k) = t_diff_R(
n,
i,
k) * t_P(
n,
j);
163 t_d_b_d_omega(
i,
j,
k) =
164 (t_d_rotated_p_d_omega(
i,
j,
k) || t_d_rotated_p_d_omega(
j,
i,
k)) / 2.;
167 t_d_u_d_omega(
i,
j,
k) = t_d_u_d_b(
i,
j,
l,
m) * t_d_b_d_omega(
l,
m,
k);
170 t_d_h_d_omega(
i,
j,
k) = t_R(
i,
l) * t_d_u_d_omega(
j,
l,
k);
172 return std::make_tuple(t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega);
198 *
dataAtPts->getStretchTensorAtPts(), nb_integration_pts);
200 *
dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
202 *
dataAtPts->getStretchH1AtPts(), nb_integration_pts);
204 *
dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
206 *
dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
208 *
dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
210 *
dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
212 *
dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
215 *
dataAtPts->getDeformationGradient(), nb_integration_pts);
217 dataAtPts->hdOmegaAtPts, nb_integration_pts);
219 dataAtPts->hdLogStretchAtPts, nb_integration_pts);
222 dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
224 dataAtPts->leviKirchhoff0AtPts, nb_integration_pts);
226 dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
228 dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
230 dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
233 dataAtPts->rotMatAtPts, nb_integration_pts);
235 *
dataAtPts->getEigenVals(), nb_integration_pts);
237 *
dataAtPts->getEigenVecs(), nb_integration_pts);
238 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
240 dataAtPts->eigenValsC, nb_integration_pts);
242 dataAtPts->eigenVecsC, nb_integration_pts);
243 dataAtPts->nbUniqC.resize(nb_integration_pts,
false);
246 dataAtPts->logStretch2H1AtPts, nb_integration_pts);
248 dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
251 dataAtPts->internalStressAtPts, nb_integration_pts);
259 auto t_levi_kirchhoff =
261 auto t_levi_kirchhoff0 =
263 auto t_levi_kirchhoff_domega =
265 auto t_levi_kirchhoff_dstreach =
267 auto t_levi_kirchhoff_dP =
269 auto t_approx_P_adjoint_dstretch =
271 auto t_approx_P_adjoint_log_du =
273 auto t_approx_P_adjoint_log_du_dP =
275 auto t_approx_P_adjoint_log_du_domega =
285 auto t_eigen_vals_C =
dataAtPts->getFTensorEigenValsC(nb_integration_pts);
286 auto t_eigen_vecs_C =
dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
293 auto t_log_stretch_total =
315 ++t_levi_kirchhoff_domega;
316 ++t_levi_kirchhoff_dstreach;
317 ++t_levi_kirchhoff_dP;
318 ++t_approx_P_adjoint_dstretch;
319 ++t_approx_P_adjoint_log_du;
320 ++t_approx_P_adjoint_log_du_dP;
321 ++t_approx_P_adjoint_log_du_domega;
334 ++t_log_stretch_total;
349 auto calculate_stretch_from_log = [&](
auto &t_log_u_src,
auto &t_u_dst,
350 auto &t_eigen_vals_dst,
351 auto &t_eigen_vecs_dst,
356 eigen_vec(
i,
j) = t_log_u_src(
i,
j);
358 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
362 nb_uniq_dst = getUniqNb<3>(eig);
363 if (nb_uniq_dst < 3) {
364 CHKERR sortEigenVals<3>(eig, eigen_vec);
366 t_eigen_vals_dst(
i) = eig(
i);
367 t_eigen_vecs_dst(
i,
j) = eigen_vec(
i,
j);
373 auto calculate_log_stretch = [&]() {
376 CHKERR calculate_stretch_from_log(t_log_u, t_u, t_eigen_vals, t_eigen_vecs,
378 t_nb_uniq = nb_uniq_val;
379 auto get_t_diff_u = [&]() {
384 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
386 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
390 auto calculate_total_stretch = [&](
auto &t_h1) {
394 t_log_u2_h1(
i,
j) = 0;
395 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
404 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
405 t_eigen_vec(
i,
j) = t_C_h1(
i,
j);
408 "Failed to compute eigenvalues of F_H1^T F_H1");
411 t_nb_uniq_C = getUniqNb<3>(t_eig_C);
412 if (t_nb_uniq_C < 3) {
413 CHKERR sortEigenVals<3>(t_eig_C, t_eigen_vec);
415 for (
int aa = 0; aa != 3; ++aa) {
416 if (!std::isfinite(t_eig_C(aa)) || t_eig_C(aa) <= 0.) {
418 "F_H1^T F_H1 must be positive definite; eigenvalue %d is "
422 const double principal_stretch = std::sqrt(t_eig_C(aa));
423 const double coordinate_stretch =
425 if (!std::isfinite(coordinate_stretch)) {
426 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
427 "Non-finite H1 coordinate stretch for principal stretch %g",
430 t_coordinate_stretch(aa) = coordinate_stretch;
432 t_eigen_vals_C(
i) = t_eig_C(
i);
433 t_eigen_vecs_C(
i,
j) = t_eigen_vec(
i,
j);
437 [](
const double v) {
return v; })(
i,
j);
440 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) + t_log_u(
i,
j);
445 auto no_h1_loop = [&]() {
455 "no_h1_loop is only implemented for LARGE_ROT");
458 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
466 CHKERR calculate_log_stretch();
469 t_u0(
i,
j) = t_u(
i,
j);
474 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
475 t_eigen_vecs_0, nb_uniq_0);
478 CHKERR calculate_total_stretch(t_h1);
480 t_u_h1(
i,
j) = t_u(
i,
j);
481 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i,
j,
k,
l);
483 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
488 auto large_rot = [&]() {
492 t_diff_diff_R(
i,
j,
k,
l) =
499 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
502 t_rotated_P(
l,
k) = t_R(
i,
l) * t_approx_P(
i,
k);
503 t_approx_P_adjoint_dstretch(
l,
k) =
504 t_diff_sym(
l,
k,
i,
j) * t_rotated_P(
i,
j);
505 t_approx_P_adjoint_log_du(
L) =
506 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
508 t_levi_kirchhoff(
m) =
509 t_diff_R(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k));
510 t_levi_kirchhoff0(
m) =
511 t_diff_R0(
i,
l,
m) * (t_u0(
l,
k) * t_approx_P0(
i,
k));
514 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
515 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
517 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
518 t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
521 t_A(
k,
l,
m) = t_diff_R(
i,
l,
m) * t_approx_P(
i,
k);
522 t_approx_P_adjoint_log_du_domega(
m,
L) =
523 t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L);
525 t_levi_kirchhoff_dstreach(
m,
L) =
526 t_diff_R(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k));
527 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
528 t_levi_kirchhoff_domega(
m,
n) =
529 t_diff_diff_R(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k));
536 auto [t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega] =
537 getDiffSpatialGradientDR(t_d_u_d_b, t_R, t_diff_R,
540 t_h_domega(
i,
k,
m) += t_d_h_d_omega(
i,
k,
m);
543 t_d_u_contract_p(
i,
l,
n) =
544 t_d_u_d_omega(
l,
k,
n) * t_approx_P(
i,
k);
545 t_levi_kirchhoff_domega(
m,
n) +=
546 t_diff_R(
i,
l,
m) * t_d_u_contract_p(
i,
l,
n);
552 t_d_b_d_p(
i,
j,
k,
l) =
553 t_diff_sym(
i,
j,
m,
l) * t_R(
k,
m);
557 t_d_u_d_p(
i,
j,
k,
l) =
558 t_d_u_d_b(
i,
j,
m,
n) * t_d_b_d_p(
m,
n,
k,
l);
559 t_levi_kirchhoff_dP(
m,
k,
l) +=
560 t_d_u_d_p(
i,
j,
k,
l) * t_d_b_d_omega(
i,
j,
m);
566 auto moderate_rot = [&](
auto &t_omega0) {
568 "moderate_rot is not implemented yet");
571 auto small_rot = [&]() {
572 t_u_h1(
i,
j) = t_u(
i,
j);
573 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i,
j,
k,
l);
575 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
577 t_R(
i,
j) =
t_kd(
i,
j) + levi_civita(
i,
j,
k) * t_omega(
k);
578 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
580 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
581 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
585 t_rotated_P(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
586 t_approx_P_adjoint_dstretch(
i,
j) =
587 t_diff_sym(
i,
j,
k,
l) * t_rotated_P(
k,
l);
588 t_approx_P_adjoint_log_du(
L) =
589 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
590 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
591 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
594 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
595 t_levi_kirchhoff0(
k) = levi_civita(
i,
j,
k) * t_approx_P0(
i,
j);
596 t_levi_kirchhoff_dstreach(
m,
L) = 0;
597 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
598 t_levi_kirchhoff_domega(
m,
n) = 0;
607 moderate_rot(t_omega0);
614 "rotationSelector not handled");
623 auto large_loop = [&]() {
633 "rotSelector should be large or small");
636 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
647 "Selected grad approximator not handled");
651 CHKERR calculate_log_stretch();
654 t_u0(
i,
j) = t_u(
i,
j);
659 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
660 t_eigen_vecs_0, nb_uniq_0);
663 CHKERR calculate_total_stretch(t_h1);
665 t_u_h1(
l,
k) = t_u(
l, o) * t_h1(o,
k);
667 t_u_h10(
l,
k) = t_u0(
l, o) * t_h1(o,
k);
668 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i, o,
k,
l) * t_h1(o,
j);
670 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
679 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
680 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
681 t_diff_R0(
i,
j,
k) = levi_civita(
i,
j,
k);
682 t_diff_diff_R(
i,
j,
l,
m) = 0;
690 t_diff_diff_R(
i,
j,
k,
l) =
696 "rotationSelector not handled");
700 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
705 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
706 t_approx_P_adjoint_dstretch(
l, o) =
707 t_diff_sym(
l, o,
i,
j) * t_rotated_P(
i,
j);
708 t_approx_P_adjoint_log_du(
L) =
709 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
712 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
713 t_levi_kirchhoff0(
m) =
714 t_diff_R0(
i,
l,
m) * t_u_h10(
l,
k) * t_approx_P0(
i,
k);
718 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
719 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
721 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
722 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
725 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
726 t_approx_P_adjoint_log_du_domega(
m,
L) =
727 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
729 t_levi_kirchhoff_dstreach(
m,
L) =
730 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
732 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
733 t_levi_kirchhoff_domega(
m,
n) =
734 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
743 auto moderate_loop = [&]() {
753 "rotSelector should be large or small");
756 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
767 "Selected grad approximator not handled");
771 CHKERR calculate_log_stretch();
773 CHKERR calculate_total_stretch(t_h1);
777 t_u_h1(
l,
k) = (
t_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
779 t_u_h10(
l,
k) = (
t_kd(
l, o) + t_log_u0(
l, o)) * t_h1(o,
k);
780 t_diff_u_h1(
i,
j,
k,
l) = t_diff(
i, o,
k,
l) * t_h1(o,
j);
782 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
791 t_R(
i,
k) =
t_kd(
i,
k) + levi_civita(
i,
k,
l) * t_omega(
l);
792 t_diff_R(
i,
j,
k) = levi_civita(
i,
j,
k);
793 t_diff_R0(
i,
j,
k) = levi_civita(
i,
j,
k);
794 t_diff_diff_R(
i,
j,
l,
m) = 0;
802 t_diff_diff_R(
i,
j,
k,
l) =
808 "rotationSelector not handled");
812 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
817 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
818 t_approx_P_adjoint_dstretch(
l, o) =
819 t_diff_sym(
l, o,
i,
j) * t_rotated_P(
i,
j);
820 t_approx_P_adjoint_log_du(
L) =
821 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
824 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
825 t_levi_kirchhoff0(
m) =
826 t_diff_R0(
i,
l,
m) * t_u_h10(
l,
k) * t_approx_P0(
i,
k);
830 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
831 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
833 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
834 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
837 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
838 t_approx_P_adjoint_log_du_domega(
m,
L) =
839 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
841 t_levi_kirchhoff_dstreach(
m,
L) =
842 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
844 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
845 t_levi_kirchhoff_domega(
m,
n) =
846 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
855 auto small_loop = [&]() {
862 "rotSelector should be small");
865 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
874 "gradApproximator not handled");
880 "stretchSelector should be linear for small loop");
882 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
883 t_u_h1(
i,
j) = t_u(
i,
j);
884 t_diff_u_h1(
i,
j,
k,
l) =
886 t_diff_u_h1(
i,
j,
k,
l) /= 2.;
887 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
889 t_log_u2_h1(
i,
j) = 0;
890 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
892 t_R(
i,
j) =
t_kd(
i,
j) + levi_civita(
i,
j,
k) * t_omega(
k);
893 t_h(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k) + t_u(
i,
j);
895 t_h_domega(
i,
j,
k) = levi_civita(
i,
j,
k);
896 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
899 t_approx_P_adjoint_dstretch(
i,
j) =
900 t_diff_sym(
i,
j,
k,
l) * t_approx_P(
k,
l);
901 t_approx_P_adjoint_log_du(
L) =
902 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
903 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
904 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
907 t_levi_kirchhoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
908 t_levi_kirchhoff0(
k) = levi_civita(
i,
j,
k) * t_approx_P0(
i,
j);
909 t_levi_kirchhoff_dstreach(
m,
L) = 0;
910 t_levi_kirchhoff_dP(
k,
i,
j) = levi_civita(
i,
j,
k);
911 t_levi_kirchhoff_domega(
m,
n) = 0;
920 case NO_H1_CONFIGURATION:
937 "gradApproximator not handled");
949 auto n_in_the_loop = getNinTheLoop();
950 auto loop_size = getLoopSize();
951 auto sense = getSkeletonSense();
952 auto nb_gauss_pts = getGaussPts().size2();
953 auto t_normal = getFTensor1NormalsAtGaussPts();
955 auto t_sigma =
dataAtPts->getFTensorApproxP(getGaussPts().size2());
959 if (!n_in_the_loop) {
963 auto t_traction = get_tracion();
964 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
966 t_sigma(
i,
j) * sense * (t_normal(
j) / t_normal.l2()) / loop_size;
984 int nb_integration_pts = getGaussPts().size2();
985 auto t_w = getFTensor0IntegrationWeight();
986 auto t_traction =
dataAtPts->getFTensorTraction(nb_integration_pts);
987 auto t_coords = getFTensor1CoordsAtGaussPts();
988 auto t_spatial_disp =
dataAtPts->getFTensorSmallWL2(nb_integration_pts);
996 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
997 double a = t_w * getMeasure();
998 loc_reaction_forces(
i) +=
a*t_traction(
i);
999 t_coords_spatial(
i) = t_coords(
i) + t_spatial_disp(
i);
1001 loc_moment_forces(
i) +=
1002 (
a * (FTensor::levi_civita<double>(
i,
j,
k) * t_coords_spatial(
j))) *
1023 int nb_integration_pts = data.
getN().size1();
1026 auto t_div_P =
dataAtPts->getFTensorDivP(nb_integration_pts);
1027 auto t_s_dot_w =
dataAtPts->getFTensorSmallWL2Dot(nb_integration_pts);
1028 auto w_l2_dot_dot_at_pts =
dataAtPts->getSmallWL2DotDotAtPts();
1029 const bool reset_w_l2_dot_dot =
1030 w_l2_dot_dot_at_pts->size1() != nb_integration_pts ||
1031 w_l2_dot_dot_at_pts->size2() != 3;
1033 *w_l2_dot_dot_at_pts, nb_integration_pts);
1034 if (reset_w_l2_dot_dot) {
1035 w_l2_dot_dot_at_pts->clear();
1037 auto t_s_dot_dot_w =
dataAtPts->getFTensorSmallWL2DotDot(nb_integration_pts);
1039 auto piola_scale =
dataAtPts->piolaScale;
1040 auto alpha_w =
alphaW / piola_scale;
1041 auto alpha_rho =
alphaRho / piola_scale;
1043 int nb_base_functions = data.
getN().size2();
1047 auto get_ftensor1 = [](
auto &
v) {
1059 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1061 auto t_nf = get_ftensor1(
nF);
1063 for (; bb != nb_dofs / 3; ++bb) {
1064 t_nf(
i) -=
a * t_row_base_fun * t_div_P(
i);
1065 t_nf(
i) +=
a * t_row_base_fun * alpha_w * t_s_dot_w(
i);
1066 t_nf(
i) +=
a * t_row_base_fun * alpha_rho * t_s_dot_dot_w(
i);
1070 for (; bb != nb_base_functions; ++bb)
1084 auto t_levi_kirchhoff =
1085 dataAtPts->getFTensorLeviKirchhoff(nb_integration_pts);
1086 auto t_levi_kirchhoff0 =
1087 dataAtPts->getFTensorLeviKirchhoff0(nb_integration_pts);
1088 auto t_omega =
dataAtPts->getFTensorRotAxis(nb_integration_pts);
1089 auto t_omega_grad =
dataAtPts->getFTensorRotAxisGrad(nb_integration_pts);
1090 auto t_omega_dot =
dataAtPts->getFTensorRotAxisDot(nb_integration_pts);
1091 auto t_omega_grad_dot =
1092 dataAtPts->getFTensorRotAxisGradDot(nb_integration_pts);
1093 int nb_base_functions = data.
getN().size2();
1099 auto get_ftensor1 = [](
auto &
v) {
1105 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1107 auto nrm_levi_kirchhoff0 = t_levi_kirchhoff0.
l2();
1110 auto t_nf = get_ftensor1(
nF);
1112 for (; bb != nb_dofs / 3; ++bb) {
1113 t_nf(
k) -= (
a * t_row_base_fun) * t_levi_kirchhoff(
k);
1116 (t_row_base_fun * t_omega(
k));
1118 (t_row_grad_fun(
i) * t_omega_grad(
k,
i));
1122 (t_row_base_fun * t_omega_dot(
k));
1126 (t_row_grad_fun(
i) * t_omega_grad_dot(
k,
i));
1131 for (; bb != nb_base_functions; ++bb) {
1141 ++t_levi_kirchhoff0;
1149 int nb_integration_pts = data.
getN().size1();
1153 int nb_base_functions = data.
getN().size2() / 3;
1161 auto get_ftensor1 = [](
auto &
v) {
1166 auto t_h =
dataAtPts->getFTensorSmallH(nb_integration_pts);
1168 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1170 auto t_nf = get_ftensor1(
nF);
1178 for (; bb != nb_dofs / 3; ++bb) {
1179 t_nf(
i) -=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
1184 for (; bb != nb_base_functions; ++bb)
1197 int nb_integration_pts = data.
getN().size1();
1201 int nb_base_functions = data.
getN().size2() / 9;
1209 auto get_ftensor0 = [](
auto &
v) {
1213 auto t_h =
dataAtPts->getFTensorSmallH(nb_integration_pts);
1215 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1217 auto t_nf = get_ftensor0(
nF);
1221 t_residuum(
i,
j) = t_h(
i,
j);
1224 for (; bb != nb_dofs; ++bb) {
1225 t_nf -=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1229 for (; bb != nb_base_functions; ++bb) {
1242 int nb_integration_pts = data.
getN().size1();
1245 auto t_w_l2 =
dataAtPts->getFTensorSmallWL2(nb_integration_pts);
1246 int nb_base_functions = data.
getN().size2() / 3;
1249 auto get_ftensor1 = [](
auto &
v) {
1254 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1256 auto t_nf = get_ftensor1(
nF);
1258 for (; bb != nb_dofs / 3; ++bb) {
1259 double div_row_base = t_row_diff_base_fun(
i,
i);
1260 t_nf(
i) -=
a * div_row_base * t_w_l2(
i);
1262 ++t_row_diff_base_fun;
1264 for (; bb != nb_base_functions; ++bb) {
1265 ++t_row_diff_base_fun;
1279 int nb_integration_pts = getGaussPts().size2();
1282 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1284 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1285 if (tag_length != 9) {
1287 "Number of internal stress components should be 9 but is %d",
1292 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1293 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(
1294 tag, &fe_ent, 1, &*const_stress_vec.data().begin());
1295 auto t_const_stress = getFTensor1FromArray<9, 9>(const_stress_vec);
1297 auto get_internal_stress =
1299 dataAtPts->internalStressAtPts, nb_integration_pts);
1300 dataAtPts->internalStressAtPts.clear();
1301 auto t_internal_stress = get_internal_stress();
1304 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1305 t_internal_stress(
L) = t_const_stress(
L);
1306 ++t_internal_stress;
1317 int nb_integration_pts = getGaussPts().size2();
1320 CHKERR getPtrFE() -> mField.get_moab().tag_get_handle(tagName.c_str(), tag);
1322 CHKERR getPtrFE() -> mField.get_moab().tag_get_length(tag, tag_length);
1323 if (tag_length != 9) {
1325 "Number of internal stress components should be 9 but is %d",
1329 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1332 CHKERR getPtrFE() -> mField.get_moab().get_connectivity(fe_ent, vert_conn,
1335 CHKERR getPtrFE() -> mField.get_moab().tag_get_data(tag, vert_conn, vert_num,
1338 auto get_internal_stress =
1340 dataAtPts->internalStressAtPts, nb_integration_pts);
1341 dataAtPts->internalStressAtPts.clear();
1342 auto t_internal_stress = get_internal_stress();
1347 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1348 auto t_vert_data = getFTensor1FromArray<9, 9>(vert_data);
1349 for (
int bb = 0; bb != nb_shape_fn; ++bb) {
1350 t_internal_stress(
L) += t_vert_data(
L) * t_shape_n;
1354 ++t_internal_stress;
1366 int nb_integration_pts = data.
getN().size1();
1367 auto v = getVolume();
1368 auto t_w = getFTensor0IntegrationWeight();
1373 auto get_ftensor2 = [](
auto &
v) {
1375 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1378 auto t_internal_stress =
1379 dataAtPts->getFTensorInternalStress(nb_integration_pts);
1383 : getFEMethod()->ts_t;
1386 double scale = scalingMethodPtr->getScale(time);
1391 int nb_base_functions = data.
getN().size2();
1393 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1395 auto t_nf = get_ftensor2(nF);
1398 t_symm_stress(
i,
j) =
1399 (t_internal_stress(
i,
j) + t_internal_stress(
j,
i)) / 2;
1402 t_residual(
L) = t_L(
i,
j,
L) * (
scale * t_symm_stress(
i,
j));
1405 for (; bb != nb_dofs / 6; ++bb) {
1406 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1410 for (; bb != nb_base_functions; ++bb)
1414 ++t_internal_stress;
1424 int nb_integration_pts = data.
getN().size1();
1425 auto v = getVolume();
1426 auto t_w = getFTensor0IntegrationWeight();
1428 auto get_ftensor2 = [](
auto &
v) {
1430 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
1433 auto t_internal_stress =
1434 dataAtPts->getFTensorInternalStressVec(nb_integration_pts);
1439 t_L = voigt_to_symm();
1443 : getFEMethod()->ts_t;
1446 double scale = scalingMethodPtr->getScale(time);
1448 int nb_base_functions = data.
getN().size2();
1450 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1452 auto t_nf = get_ftensor2(nF);
1455 t_residual(
L) = t_L(
M,
L) * (
scale * t_internal_stress(
M));
1458 for (; bb != nb_dofs / 6; ++bb) {
1459 t_nf(
L) +=
a * t_row_base_fun * t_residual(
L);
1463 for (; bb != nb_base_functions; ++bb)
1467 ++t_internal_stress;
1472template <AssemblyType A>
1478 for (
auto &bc : (*bcDispPtr)) {
1480 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1483 int nb_integration_pts = OP::getGaussPts().size2();
1484 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1485 auto t_w = OP::getFTensor0IntegrationWeight();
1486 int nb_base_functions = data.
getN().size2() / 3;
1493 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1495 scale *= scalingMethodsMap.at(bc.blockName)
1498 scale *= scalingMethodsMap.at(bc.blockName)
1499 ->getScale(OP::getFEMethod()->ts_t);
1503 <<
"No scaling method found for " << bc.blockName;
1510 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1511 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1513 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1515 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
1519 for (; bb != nb_base_functions; ++bb)
1531 return OP::iNtegrate(data);
1539 for (
auto &bc : (*bcDispPtr)) {
1541 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1544 int nb_integration_pts = OP::getGaussPts().size2();
1545 auto t_w = OP::getFTensor0IntegrationWeight();
1546 int nb_base_functions = data.
getN().size2();
1549 if (!this->sourceVec) {
1551 "Source vector for OpTauStabilizationDispRhsBc is not set");
1553 if (data.
getN().size1() != nb_integration_pts) {
1555 "Number of integration points in data should be %d but is %d",
1556 nb_integration_pts, (
int)data.
getN().size1());
1558 if (nb_base_functions < nb_dofs /
SPACE_DIM) {
1560 "Number of base functions in data should be %d but is %d",
1568 *this->sourceVec, nb_integration_pts)();
1580 ->getScale(OP::getFEMethod()->ts_t);
1584 <<
"No scaling method found for " << bc.blockName;
1591 auto area = getMeasure();
1592 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1593 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1595 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1599 bc.flags[ii] ? t_disp_val(ii) - t_bc_disp(ii) : 0.;
1600 auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
1602 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1604 (tau_scale * t_row_base_fun) * t_bc_residual(
i);
1608 for (; bb != nb_base_functions; ++bb)
1627 for (
auto &bc : (*bcDispPtr)) {
1629 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1632 int nb_integration_pts = OP::getGaussPts().size2();
1633 auto t_w = OP::getFTensor0IntegrationWeight();
1634 int nb_base_functions = row_data.
getN().size2();
1642 t_bc_mask(ii) = bc.flags[ii] ? 1. : 0.;
1644 auto get_t_vec = [&](
const int rr) {
1645 std::array<double *, SPACE_DIM> ptrs;
1647 ptrs[
i] = &OP::locMat(rr +
i,
i);
1652 auto area = getMeasure();
1653 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1654 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1656 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1658 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1661 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
1663 tau_scale * (t_row_base_fun * t_col_base_fun) * t_bc_mask(
i);
1669 for (; rr != nb_base_functions; ++rr)
1685 for (
auto &bc : (*bcDispPtr)) {
1686 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1688 auto analytical_data = getAnalyticalExpr(
this,
analytical_expr, bc.blockName);
1689 auto &v_analytical_expr = std::get<1>(analytical_data);
1692 int nb_integration_pts = OP::getGaussPts().size2();
1693 auto t_w = OP::getFTensor0IntegrationWeight();
1694 int nb_base_functions = data.
getN().size2();
1698 if (!this->sourceVec) {
1700 "Source vector for OpTauStabilizationOpAnalyticalDispBc is not "
1703 if (data.
getN().size1() != nb_integration_pts) {
1705 "Number of integration points in data should be %d but is %d",
1706 nb_integration_pts, (
int)data.
getN().size1());
1708 if (nb_base_functions < nb_dofs /
SPACE_DIM) {
1710 "Number of base functions in data should be at least %d but is "
1712 nb_dofs /
SPACE_DIM, nb_base_functions);
1718 *this->sourceVec, nb_integration_pts)();
1723 auto area = getMeasure();
1724 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1725 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1727 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1731 bc.flags[ii] ? t_disp_val(ii) - t_bc_disp(ii) : 0.;
1732 auto t_nf = getFTensor1FromPtr<3, 3>(OP::locF.data().data());
1734 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1736 (tau_scale * t_row_base_fun) * t_bc_residual(
i);
1740 for (; bb != nb_base_functions; ++bb)
1759 for (
auto &bc : (*bcDispPtr)) {
1760 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1762 int nb_integration_pts = OP::getGaussPts().size2();
1763 auto t_w = OP::getFTensor0IntegrationWeight();
1764 int nb_base_functions = row_data.
getN().size2();
1770 t_bc_mask(ii) = bc.flags[ii] ? 1. : 0.;
1772 auto get_t_vec = [&](
const int rr) {
1773 std::array<double *, SPACE_DIM> ptrs;
1775 ptrs[
i] = &OP::locMat(rr +
i,
i);
1780 auto area = getMeasure();
1781 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1782 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1784 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1786 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
1789 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
1791 tau_scale * (t_row_base_fun * t_col_base_fun) * t_bc_mask(
i);
1797 for (; rr != nb_base_functions; ++rr)
1809template <AssemblyType A>
1817 double time = OP::getFEMethod()->ts_t;
1825 for (
auto &bc : (*bcRotPtr)) {
1827 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1829 int nb_integration_pts = OP::getGaussPts().size2();
1830 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
1831 auto t_w = OP::getFTensor0IntegrationWeight();
1833 int nb_base_functions = data.
getN().size2() / 3;
1836 auto get_ftensor1 = [](
auto &
v) {
1849 auto get_rotation_angle = [&]() {
1850 double theta = bc.theta;
1851 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
1852 theta *= scalingMethodsMap.at(bc.blockName)->getScale(time);
1857 auto get_rotation = [&](
auto theta) {
1859 if (bc.vals.size() == 7) {
1860 t_omega(0) = bc.vals[4];
1861 t_omega(1) = bc.vals[5];
1862 t_omega(2) = bc.vals[6];
1865 t_omega(
i) = OP::getFTensor1Normal()(
i);
1867 if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
1871 <<
"Rotation axis is zero vector for block " << bc.blockName
1872 <<
". This may lead to unexpected results.";
1874 t_omega(
i) *= theta;
1876 RotSelector::SMALL_ROT
1881 auto t_R = get_rotation(get_rotation_angle());
1882 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1884 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1886 t_delta(
i) = t_center(
i) - t_coords(
i);
1888 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1890 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1892 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1893 t_nf(
i) += t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1897 for (; bb != nb_base_functions; ++bb)
1910 return OP::iNtegrate(data);
1919 double time = OP::getFEMethod()->ts_t;
1927 for (
auto &bc : (*bcRotPtr)) {
1929 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1931 int nb_integration_pts = OP::getGaussPts().size2();
1932 auto t_w = OP::getFTensor0IntegrationWeight();
1934 int nb_base_functions = data.
getN().size2();
1937 auto get_ftensor1 = [](
auto &
v) {
1950 auto get_rotation_angle = [&]() {
1951 double theta = bc.theta;
1958 auto get_rotation = [&](
auto theta) {
1960 if (bc.vals.size() == 7) {
1961 t_omega(0) = bc.vals[4];
1962 t_omega(1) = bc.vals[5];
1963 t_omega(2) = bc.vals[6];
1966 t_omega(
i) = OP::getFTensor1Normal()(
i);
1968 if (t_omega.
l2() > std::numeric_limits<double>::epsilon()) {
1971 t_omega(
i) *= theta;
1973 RotSelector::SMALL_ROT
1978 auto area = getMeasure();
1979 auto t_R = get_rotation(get_rotation_angle());
1980 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
1983 *this->sourceVec, nb_integration_pts)();
1985 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1987 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1990 t_delta(
i) = t_center(
i) - t_coords(
i);
1992 t_bc_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1994 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
1996 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1998 (tau_scale * t_row_base_fun) * (t_disp_val(
i) - t_bc_disp(
i));
2002 for (; bb != nb_base_functions; ++bb)
2021 for (
auto &bc : (*bcRotPtr)) {
2023 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2026 int nb_integration_pts = OP::getGaussPts().size2();
2027 auto t_w = OP::getFTensor0IntegrationWeight();
2028 int nb_base_functions = row_data.
getN().size2();
2034 auto get_t_vec = [&](
const int rr) {
2035 std::array<double *, SPACE_DIM> ptrs;
2037 ptrs[
i] = &OP::locMat(rr +
i,
i);
2042 auto area = getMeasure();
2043 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2044 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2046 area * t_w * OP::betaCoeff(t_coords(0), t_coords(1), t_coords(2));
2048 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2051 for (
int cc = 0; cc != nb_dofs /
SPACE_DIM; ++cc) {
2052 for (
int ii = 0; ii !=
SPACE_DIM; ++ii) {
2053 t_mat(ii) += tau_scale * (t_row_base_fun * t_col_base_fun);
2060 for (; rr != nb_base_functions; ++rr)
2072template <AssemblyType A>
2076 double time = OP::getFEMethod()->ts_t;
2084 for (
auto &bc : (*bcDispPtr)) {
2086 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2088 for (
auto &bd : (*brokenBaseSideDataPtr)) {
2092 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2093 auto t_w = OP::getFTensor0IntegrationWeight();
2101 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2102 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2105 <<
"No scaling method found for " << bc.blockName;
2109 double val =
scale * bc.val;
2112 int nb_integration_pts = OP::getGaussPts().size2();
2113 int nb_base_functions = data.
getN().size2();
2115 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2118 t_N(
i) = t_normal(
i);
2122 t_P(
i,
j) = t_N(
i) * t_N(
j);
2127 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
2131 t_Q(
i,
j) * t_traction(
j) + t_P(
i,
j) * 2 * t_u(
j) - t_N(
i) * val;
2133 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
2135 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
2136 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
2140 for (; bb != nb_base_functions; ++bb)
2154template <AssemblyType A>
2160 double time = OP::getFEMethod()->ts_t;
2165 int row_nb_dofs = row_data.
getIndices().size();
2166 int col_nb_dofs = col_data.
getIndices().size();
2167 auto &locMat = OP::locMat;
2168 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2174 for (
auto &bc : (*bcDispPtr)) {
2176 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2178 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2179 auto t_w = OP::getFTensor0IntegrationWeight();
2185 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2186 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2189 <<
"No scaling method found for " << bc.blockName;
2192 int nb_integration_pts = OP::getGaussPts().size2();
2193 int row_nb_dofs = row_data.
getIndices().size();
2194 int col_nb_dofs = col_data.
getIndices().size();
2195 int nb_base_functions = row_data.
getN().size2();
2200 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2203 t_N(
i) = t_normal(
i);
2207 t_P(
i,
j) = t_N(
i) * t_N(
j);
2210 t_d_res(
i,
j) = 2.0 * t_P(
i,
j);
2213 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2214 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2217 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2218 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
2225 for (; rr != nb_base_functions; ++rr)
2232 locMat *= OP::getMeasure();
2238template <AssemblyType A>
2244 double time = OP::getFEMethod()->ts_t;
2249 int row_nb_dofs = row_data.
getIndices().size();
2250 int col_nb_dofs = col_data.
getIndices().size();
2251 auto &locMat = OP::locMat;
2252 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2258 for (
auto &bc : (*bcDispPtr)) {
2260 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2262 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2263 auto t_w = OP::getFTensor0IntegrationWeight();
2272 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2273 scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2276 <<
"No scaling method found for " << bc.blockName;
2279 int nb_integration_pts = OP::getGaussPts().size2();
2280 int nb_base_functions = row_data.
getN().size2();
2283 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2286 t_N(
i) = t_normal(
i);
2290 t_P(
i,
j) = t_N(
i) * t_N(
j);
2295 t_d_res(
i,
j) = t_Q(
i,
j);
2298 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2299 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2302 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2304 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * t_d_res(
i,
j);
2311 for (; rr != nb_base_functions; ++rr)
2318 locMat *= OP::getMeasure();
2325 return OP::iNtegrate(data);
2344 for (
auto &bc : (*bcSpringPtr)) {
2346 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2348 for (
auto &bd : (*brokenBaseSideDataPtr)) {
2352 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2353 auto t_w = OP::getFTensor0IntegrationWeight();
2361 int nb_integration_pts = OP::getGaussPts().size2();
2362 int nb_base_functions = data.
getN().size2();
2364 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2367 t_N(
i) = t_normal(
i);
2371 t_P(
i,
j) = t_N(
i) * t_N(
j);
2376 t_traction(
i) = t_approx_P(
i,
j) * t_N(
j);
2379 t_res(
i) = 0.5 *(t_traction(
i)) - bc.normalStiffness * t_P(
i,
j) * t_u(
j) -
2380 bc.tangentialStiffness * t_Q(
i,
j) * t_u(
j);
2382 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
2384 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
2385 t_nf(
i) += (t_w * t_row_base * OP::getMeasure()) * t_res(
i);
2389 for (; bb != nb_base_functions; ++bb)
2407 int row_nb_dofs = row_data.
getIndices().size();
2408 int col_nb_dofs = col_data.
getIndices().size();
2409 auto &locMat = OP::locMat;
2410 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2416 for (
auto &bc : (*bcSpringPtr)) {
2418 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2420 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2421 auto t_w = OP::getFTensor0IntegrationWeight();
2426 int nb_integration_pts = OP::getGaussPts().size2();
2427 int nb_base_functions = row_data.
getN().size2();
2432 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2435 t_N(
i) = t_normal(
i);
2439 t_P(
i,
j) = t_N(
i) * t_N(
j);
2444 t_d_res(
i,
j) = -(bc.normalStiffness * t_P(
i,
j) +
2445 bc.tangentialStiffness * t_Q(
i,
j));
2448 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2449 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2452 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2453 t_mat(
i,
j) += (t_w * t_row_base * t_col_base) * t_d_res(
i,
j);
2460 for (; rr != nb_base_functions; ++rr)
2467 locMat *= OP::getMeasure();
2477 int row_nb_dofs = row_data.
getIndices().size();
2478 int col_nb_dofs = col_data.
getIndices().size();
2479 auto &locMat = OP::locMat;
2480 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2486 for (
auto &bc : (*bcSpringPtr)) {
2488 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2490 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2491 auto t_w = OP::getFTensor0IntegrationWeight();
2497 int nb_integration_pts = OP::getGaussPts().size2();
2498 int nb_base_functions = row_data.
getN().size2();
2503 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2506 t_N(
i) = t_normal(
i);
2510 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2511 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2514 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2516 ((t_w * t_row_base) * (t_N(
k) * t_col_base(
k))) * 0.5 *
t_kd(
i,
j);
2523 for (; rr != nb_base_functions; ++rr)
2530 locMat *= OP::getMeasure();
2536template <AssemblyType A>
2540 double time = OP::getFEMethod()->ts_t;
2548 for (
auto &bc : (*bcDispPtr)) {
2550 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2555 auto [block_name, v_analytical_expr] =
2563 int nb_integration_pts = OP::getGaussPts().size2();
2564 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
2565 auto t_w = OP::getFTensor0IntegrationWeight();
2566 int nb_base_functions = data.
getN().size2() / 3;
2575 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2576 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
2579 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
2581 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
2585 for (; bb != nb_base_functions; ++bb)
2598 return OP::iNtegrate(data);
2607 int nb_integration_pts = getGaussPts().size2();
2608 int nb_base_functions = data.
getN().size2();
2610 double time = getFEMethod()->ts_t;
2616 if (this->locF.size() != nb_dofs)
2618 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2621 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2624 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
2626 auto t_w = getFTensor0IntegrationWeight();
2627 auto t_coords = getFTensor1CoordsAtGaussPts();
2628 auto t_normal = getFTensor1NormalsAtGaussPts();
2632 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2634 double a = sqrt(t_normal(
i) * t_normal(
i));
2636 const auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2637 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2639 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2641 (time_scale *
a * t_w * t_row_base * tau) * (t_val(
i) *
scale);
2646 for (; rr != nb_base_functions; ++rr)
2657 for (
auto &bc : *(
bcData)) {
2658 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2660 double time_scale = 1;
2668 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*"))) {
2672 return -y * (y - 1) / 0.25;
2674 CHKERR integrate_rhs(bc, calc_tau, time_scale);
2677 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2691 int nb_integration_pts = getGaussPts().size2();
2692 int nb_base_functions = data.
getN().size2();
2694 double time = getFEMethod()->ts_t;
2700 if (this->locF.size() != nb_dofs)
2702 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2705 auto integrate_rhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2710 auto t_w = getFTensor0IntegrationWeight();
2711 auto t_coords = getFTensor1CoordsAtGaussPts();
2712 auto t_tangent1 = getFTensor1Tangent1AtGaussPts();
2713 auto t_tangent2 = getFTensor1Tangent2AtGaussPts();
2719 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2731 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) * t_tangent1(
j)) *
2734 t_normal(
i) = (FTensor::levi_civita<double>(
i,
j,
k) *
2735 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2736 (t_tangent2(
k) + t_grad_gamma_u(
k, N1));
2738 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2740 t_val(
i) = (time_scale * t_w * tau *
scale * val) * t_normal(
i);
2742 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2744 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2745 t_f(
i) += t_row_base * t_val(
i);
2750 for (; rr != nb_base_functions; ++rr)
2765 for (
auto &bc : *(
bcData)) {
2766 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2768 double time_scale = 1;
2776 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2783template <AssemblyType A>
2794 double time = OP::getFEMethod()->ts_t;
2799 int nb_base_functions = row_data.
getN().size2();
2800 int row_nb_dofs = row_data.
getIndices().size();
2801 int col_nb_dofs = col_data.
getIndices().size();
2802 int nb_integration_pts = OP::getGaussPts().size2();
2803 auto &locMat = OP::locMat;
2804 locMat.resize(row_nb_dofs, col_nb_dofs,
false);
2807 auto integrate_lhs = [&](
auto &bc,
auto calc_tau,
double time_scale) {
2812 auto t_w = OP::getFTensor0IntegrationWeight();
2813 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
2814 auto t_tangent1 = OP::getFTensor1Tangent1AtGaussPts();
2815 auto t_tangent2 = OP::getFTensor1Tangent2AtGaussPts();
2820 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2830 auto tau = calc_tau(t_coords(0), t_coords(1), t_coords(2));
2831 auto t_val = time_scale * t_w * tau * val;
2834 for (; rr != row_nb_dofs /
SPACE_DIM; ++rr) {
2835 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
2838 for (
auto cc = 0; cc != col_nb_dofs /
SPACE_DIM; ++cc) {
2840 t_normal_du(
i,
l) = (FTensor::levi_civita<double>(
i,
j,
k) *
2841 (t_tangent2(
k) + t_grad_gamma_u(
k, N1))) *
2842 t_kd(
j,
l) * t_diff_col_base(N0)
2846 (FTensor::levi_civita<double>(
i,
j,
k) *
2847 (t_tangent1(
j) + t_grad_gamma_u(
j, N0))) *
2848 t_kd(
k,
l) * t_diff_col_base(N1);
2850 t_mat(
i,
j) += t_row_base * t_val * t_normal_du(
i,
j);
2857 for (; rr != nb_base_functions; ++rr)
2873 for (
auto &bc : *(bcData)) {
2874 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2876 double time_scale = 1;
2877 if (scalingMethodsMap.find(bc.blockName) != scalingMethodsMap.end()) {
2878 time_scale *= scalingMethodsMap.at(bc.blockName)->getScale(time);
2884 bc, [](
double,
double,
double) {
return 1; }, time_scale);
2903 int nb_integration_pts = getGaussPts().size2();
2904 int nb_base_functions = data.
getN().size2();
2906 double time = getFEMethod()->ts_t;
2912 if (this->locF.size() != nb_dofs)
2914 "Size of locF %ld != nb_dofs %d", this->locF.size(), nb_dofs);
2919 for (
auto &bc : *(
bcData)) {
2920 if (bc.faces.find(fe_ent) != bc.faces.end()) {
2924 auto [block_name, v_analytical_expr] =
2928 auto t_w = getFTensor0IntegrationWeight();
2929 auto t_coords = getFTensor1CoordsAtGaussPts();
2933 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2935 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
2937 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
2938 t_f(
i) -= t_w * t_row_base * (t_val(
i) *
scale);
2943 for (; rr != nb_base_functions; ++rr)
2949 this->locF *= getMeasure();
2958 int nb_integration_pts = row_data.
getN().size1();
2959 int row_nb_dofs = row_data.
getIndices().size();
2960 int col_nb_dofs = col_data.
getIndices().size();
2961 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
2963 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
2968 int row_nb_base_functions = row_data.
getN().size2();
2970 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2973 for (; rr != row_nb_dofs / 3; ++rr) {
2975 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2976 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2977 double div_col_base = t_col_diff_base_fun(
i,
i);
2978 t_m(
i) -=
a * t_row_base_fun * div_col_base;
2980 ++t_col_diff_base_fun;
2984 for (; rr != row_nb_base_functions; ++rr)
2995 if (
alphaW < std::numeric_limits<double>::epsilon() &&
2996 alphaRho < std::numeric_limits<double>::epsilon())
2999 const int nb_integration_pts = row_data.
getN().size1();
3000 const int row_nb_dofs = row_data.
getIndices().size();
3001 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3003 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
3012 auto piola_scale =
dataAtPts->piolaScale;
3013 auto alpha_w =
alphaW / piola_scale;
3014 auto alpha_rho =
alphaRho / piola_scale;
3016 int row_nb_base_functions = row_data.
getN().size2();
3019 double ts_scale = alpha_w *
getTSa();
3020 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon())
3021 ts_scale += alpha_rho *
getTSaa();
3023 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3024 double a =
v * t_w * ts_scale;
3027 for (; rr != row_nb_dofs / 3; ++rr) {
3030 auto t_m = get_ftensor1(
K, 3 * rr, 0);
3031 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
3032 const double b =
a * t_row_base_fun * t_col_base_fun;
3041 for (; rr != row_nb_base_functions; ++rr)
3060 int nb_integration_pts = row_data.
getN().size1();
3061 int row_nb_dofs = row_data.
getIndices().size();
3062 int col_nb_dofs = col_data.
getIndices().size();
3063 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
3066 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3068 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3070 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
3072 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
3074 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
3076 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
3082 int row_nb_base_functions = row_data.
getN().size2();
3085 auto t_approx_P_adjoint_log_du_dP =
3086 dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
3088 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3091 for (; rr != row_nb_dofs / 6; ++rr) {
3094 auto t_m = get_ftensor3(
K, 6 * rr, 0);
3096 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3098 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
3106 for (; rr != row_nb_base_functions; ++rr)
3109 ++t_approx_P_adjoint_log_du_dP;
3125 int nb_integration_pts = row_data.
getN().size1();
3126 int row_nb_dofs = row_data.
getIndices().size();
3127 int col_nb_dofs = col_data.
getIndices().size();
3128 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3130 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
3138 int row_nb_base_functions = row_data.
getN().size2();
3140 auto t_approx_P_adjoint_log_du_dP =
3141 dataAtPts->getFTensorAdjointPdUdP(nb_integration_pts);
3143 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3146 for (; rr != row_nb_dofs / 6; ++rr) {
3147 auto t_m = get_ftensor2(
K, 6 * rr, 0);
3148 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
3149 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
3151 a * (t_approx_P_adjoint_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
3158 for (; rr != row_nb_base_functions; ++rr)
3161 ++t_approx_P_adjoint_log_du_dP;
3174 int row_nb_dofs = row_data.
getIndices().size();
3175 int col_nb_dofs = col_data.
getIndices().size();
3176 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
3179 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3181 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3183 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
3185 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
3187 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
3189 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
3201 auto t_approx_P_adjoint_log_du_domega =
3202 dataAtPts->getFTensorAdjointPdUdOmega(nb_integration_pts);
3204 int row_nb_base_functions = row_data.
getN().size2();
3207 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3211 for (; rr != row_nb_dofs / 6; ++rr) {
3213 auto t_m = get_ftensor3(
K, 6 * rr, 0);
3214 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3215 double v =
a * t_row_base_fun * t_col_base_fun;
3216 t_m(
L,
k) -=
v * t_approx_P_adjoint_log_du_domega(
k,
L);
3223 for (; rr != row_nb_base_functions; ++rr)
3227 ++t_approx_P_adjoint_log_du_domega;
3237 int row_nb_dofs = row_data.
getIndices().size();
3238 int col_nb_dofs = col_data.
getIndices().size();
3239 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3243 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
3244 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
3246 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
3247 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
3249 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
3250 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5)
3260 auto t_levi_kirchhoff_du =
3261 dataAtPts->getFTensorLeviKirchhoffdLogStretch(nb_integration_pts);
3262 int row_nb_base_functions = row_data.
getN().size2();
3264 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3267 for (; rr != row_nb_dofs / 3; ++rr) {
3268 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3269 const double b =
a * t_row_base_fun;
3271 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
3272 t_m(
k,
L) -= (b * t_col_base_fun) * t_levi_kirchhoff_du(
k,
L);
3278 for (; rr != row_nb_base_functions; ++rr) {
3282 ++t_levi_kirchhoff_du;
3299 int row_nb_dofs = row_data.
getIndices().size();
3300 int col_nb_dofs = col_data.
getIndices().size();
3301 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3304 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3306 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3308 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3316 int row_nb_base_functions = row_data.
getN().size2();
3318 auto t_levi_kirchhoff_dP =
3319 dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
3321 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3324 for (; rr != row_nb_dofs / 3; ++rr) {
3325 double b =
a * t_row_base_fun;
3327 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3328 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3329 t_m(
m,
i) -= b * (t_levi_kirchhoff_dP(
m,
i,
k) * t_col_base_fun(
k));
3335 for (; rr != row_nb_base_functions; ++rr) {
3340 ++t_levi_kirchhoff_dP;
3349 int row_nb_dofs = row_data.
getIndices().size();
3350 int col_nb_dofs = col_data.
getIndices().size();
3352 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3354 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
3363 auto t_levi_kirchoff_dP =
3364 dataAtPts->getFTensorLeviKirchhoffP(nb_integration_pts);
3366 int row_nb_base_functions = row_data.
getN().size2();
3369 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3372 for (; rr != row_nb_dofs / 3; ++rr) {
3373 double b =
a * t_row_base_fun;
3374 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
3375 auto t_m = get_ftensor1(
K, 3 * rr, 0);
3376 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
3377 t_m(
m) -= b * (t_levi_kirchoff_dP(
m,
i,
k) * t_col_base_fun(
i,
k));
3384 for (; rr != row_nb_base_functions; ++rr) {
3388 ++t_levi_kirchoff_dP;
3397 int row_nb_dofs = row_data.
getIndices().size();
3398 int col_nb_dofs = col_data.
getIndices().size();
3399 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3402 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3404 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3406 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3421 (std::abs(
alphaViscousR) > std::numeric_limits<double>::epsilon() ||
3422 std::abs(
alphaViscousR0) > std::numeric_limits<double>::epsilon() ||
3428 auto t_levi_kirchhoff0 =
3429 dataAtPts->getFTensorLeviKirchhoff0(nb_integration_pts);
3430 auto t_levi_kirchhoff_domega =
3431 dataAtPts->getFTensorLeviKirchhoffdOmega(nb_integration_pts);
3432 int row_nb_base_functions = row_data.
getN().size2();
3437 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3438 auto nrm_levi_kirchhoff0 = t_levi_kirchhoff0.l2();
3450 for (; rr != row_nb_dofs / 3; ++rr) {
3451 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3452 const double row_mass =
a * t_row_base_fun;
3455 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3456 t_m(
k,
l) -= (row_mass * t_col_base_fun) * t_levi_kirchhoff_domega(
k,
l);
3457 t_m(
k,
l) +=
t_kd(
k,
l) * (mass_coeff * t_row_base_fun * t_col_base_fun);
3459 t_kd(
k,
l) * (grad_coeff * (t_row_grad_fun(
i) * t_col_grad_fun(
i)));
3467 for (; rr != row_nb_base_functions; ++rr) {
3472 ++t_levi_kirchhoff0;
3473 ++t_levi_kirchhoff_domega;
3478template <
typename TInvD,
typename TRotation>
3485 t_d_b_d_p(
i,
j,
k,
l) = t_diff_sym(
i,
j,
m,
l) * t_R(
k,
m);
3489 t_d_u_d_p(
i,
j,
k,
l) =
3490 t_d_u_d_b(
i,
j,
m,
n) * t_d_b_d_p(
m,
n,
k,
l);
3494 t_d_h_d_p(
i,
j,
k,
l) = t_R(
i,
m) * t_d_u_d_p(
m,
j,
k,
l);
3505 integrateImpl<size_symm * size_symm>(row_data, col_data));
3517 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3520 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3522 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3524 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3530 int row_nb_dofs = row_data.
getIndices().size();
3531 int col_nb_dofs = col_data.
getIndices().size();
3535 int row_nb_base_functions = row_data.
getN().size2() / 3;
3543 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(
dataAtPts->matInvD);
3544 auto t_R =
dataAtPts->getFTensorRotMat(nb_integration_pts);
3547 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3550 auto assemble = [&](
auto &t_diff_h_p) {
3552 for (; rr != row_nb_dofs / 3; ++rr) {
3554 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3555 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3557 a * t_row_base(
j) * (t_diff_h_p(
i,
j,
k,
l) * t_col_base(
l));
3565 for (; rr != row_nb_base_functions; ++rr)
3571 auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
3592 integrateImpl<size_symm * size_symm>(row_data, col_data));
3606 int row_nb_dofs = row_data.
getIndices().size();
3607 int col_nb_dofs = col_data.
getIndices().size();
3611 int row_nb_base_functions = row_data.
getN().size2() / 9;
3619 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(
dataAtPts->matInvD);
3620 auto t_R =
dataAtPts->getFTensorRotMat(nb_integration_pts);
3623 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3626 auto assemble = [&](
auto &t_diff_h_p) {
3628 for (; rr != row_nb_dofs; ++rr) {
3630 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
3631 K(rr, cc) -=
a * (t_row_base(
i,
j) *
3632 (t_diff_h_p(
i,
j,
k,
l) * t_col_base(
k,
l)));
3639 for (; rr != row_nb_base_functions; ++rr)
3645 auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
3665 integrateImpl<size_symm * size_symm>(row_data, col_data));
3678 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3681 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
3687 int row_nb_dofs = row_data.
getIndices().size();
3688 int col_nb_dofs = col_data.
getIndices().size();
3692 int row_nb_base_functions = row_data.
getN().size2() / 9;
3702 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(
dataAtPts->matInvD);
3703 auto t_R =
dataAtPts->getFTensorRotMat(nb_integration_pts);
3706 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3709 auto assemble = [&](
auto &t_diff_h_p) {
3710 auto t_m = get_ftensor1(
K, 0, 0);
3712 for (; rr != row_nb_dofs; ++rr) {
3714 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3716 a * (t_row_base(
i,
j) * t_diff_h_p(
i,
j,
k,
l)) *
3725 for (; rr != row_nb_base_functions; ++rr)
3731 auto t_diff_h_p = getDiffSpatialGradientDP(t_inv_D, t_R);
3748 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
3751 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2)
3757 int row_nb_dofs = row_data.
getIndices().size();
3758 int col_nb_dofs = col_data.
getIndices().size();
3762 int row_nb_base_functions = row_data.
getN().size2() / 9;
3768 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3771 auto t_m = get_ftensor1(
K, 0, 0);
3774 for (; rr != row_nb_dofs; ++rr) {
3776 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3777 t_m(
k) +=
a * t_row_base(
k,
l) * t_col_base(
l);
3785 for (; rr != row_nb_base_functions; ++rr)
3804 int nb_integration_pts = row_data.
getN().size1();
3805 int row_nb_dofs = row_data.
getIndices().size();
3806 int col_nb_dofs = col_data.
getIndices().size();
3808 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3811 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
3813 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
3815 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
3822 int row_nb_base_functions = row_data.
getN().size2() / 3;
3825 auto t_h_domega =
dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
3827 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3831 for (; rr != row_nb_dofs / 3; ++rr) {
3834 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
3837 auto t_m = get_ftensor2(
K, 3 * rr, 0);
3838 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3839 t_m(
i,
j) -= (
a * t_col_base_fun) * t_PRT(
i,
j);
3847 for (; rr != row_nb_base_functions; ++rr)
3867 int nb_integration_pts = row_data.
getN().size1();
3868 int row_nb_dofs = row_data.
getIndices().size();
3869 int col_nb_dofs = col_data.
getIndices().size();
3871 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
3873 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
3878 int row_nb_base_functions = row_data.
getN().size2() / 9;
3881 auto t_h_domega =
dataAtPts->getFTensorSmallHdOmega(nb_integration_pts);
3882 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
3886 for (; rr != row_nb_dofs; ++rr) {
3889 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
3892 auto t_m = get_ftensor2(
K, rr, 0);
3893 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
3894 t_m(
j) -= (
a * t_col_base_fun) * t_PRT(
j);
3902 for (; rr != row_nb_base_functions; ++rr)
3915 if (
tagSense != getSkeletonSense())
3918 auto get_tag = [&](
auto name) {
3919 auto &mob = getPtrFE()->mField.get_moab();
3925 auto get_tag_value = [&](
auto &&tag,
int dim) {
3926 auto &mob = getPtrFE()->mField.get_moab();
3927 auto face = getSidePtrFE()->getFEEntityHandle();
3928 std::vector<double> value(dim);
3929 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3933 auto create_tag = [
this](
const std::string tag_name,
const int size) {
3934 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3937 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3942 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3943 Tag th_detF = create_tag(
"detF", 1);
3944 Tag th_traction = create_tag(
"traction", 3);
3945 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3947 Tag th_energy = create_tag(
"Energy", 1);
3948 Tag th_young_modulus = create_tag(
"YoungModulus", 1);
3950 const auto nb_gauss_pts = getGaussPts().size2();
3951 auto t_w =
dataAtPts->getFTensorSmallWL2(nb_gauss_pts);
3952 auto t_h =
dataAtPts->getFTensorSmallH(nb_gauss_pts);
3953 auto t_approx_P =
dataAtPts->getFTensorApproxP(nb_gauss_pts);
3955 auto t_normal = getFTensor1NormalsAtGaussPts();
3956 auto t_disp =
dataAtPts->getFTensorSmallWH1(nb_gauss_pts);
3960 if (
dataAtPts->energyAtPts.size() == 0) {
3962 dataAtPts->energyAtPts.resize(nb_gauss_pts);
3983 auto set_float_precision = [](
const double x) {
3984 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3991 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3993 v = set_float_precision(
v);
4001 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
4004 for (
auto &
a :
v.data())
4005 a = set_float_precision(
a);
4007 &*
v.data().begin());
4015 &
m(0, 0), &
m(0, 1), &
m(0, 2),
4017 &
m(1, 0), &
m(1, 1), &
m(1, 2),
4019 &
m(2, 0), &
m(2, 1), &
m(2, 2));
4021 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
4023 t_m(
i,
j) = t_d(
i,
j);
4024 for (
auto &
v :
m.data())
4025 v = set_float_precision(
v);
4027 &*
m.data().begin());
4031 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
4034 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
4037 CHKERR save_vec_tag(th_traction, t_traction, gg);
4039 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
4040 if (!std::isfinite(u_error))
4042 CHKERR save_scal_tag(th_disp_error, u_error, gg);
4043 CHKERR save_scal_tag(th_energy, t_energy, gg);
4045 CHKERR save_scal_tag(th_young_modulus, t_youngs_modulus, gg);
4049 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
4050 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
4060 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
4061 std::vector<FieldSpace> spaces, std::string geom_field_name,
4062 boost::shared_ptr<Range> crack_front_edges_ptr) {
4065 constexpr bool scale_l2 =
false;
4069 "Scale L2 Ainsworth Legendre base is not implemented");
4078 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
4079 std::vector<FieldSpace> spaces, std::string geom_field_name,
4080 boost::shared_ptr<Range> crack_front_edges_ptr) {
4083 constexpr bool scale_l2 =
false;
4087 "Scale L2 Ainsworth Legendre base is not implemented");
4096 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
4097 std::vector<FieldSpace> spaces, std::string geom_field_name,
4098 boost::shared_ptr<Range> crack_front_edges_ptr,
4099 boost::shared_ptr<MatrixDouble> jac, boost::shared_ptr<VectorDouble> det,
4100 boost::shared_ptr<MatrixDouble> inv_jac) {
4103 if (!geom_field_name.empty()) {
4105 auto jac = boost::make_shared<MatrixDouble>();
4106 auto det = boost::make_shared<VectorDouble>();
4109 geom_field_name, jac));
4116 constexpr bool scale_l2_ainsworth_legendre_base =
false;
4118 if (scale_l2_ainsworth_legendre_base) {
4126 boost::shared_ptr<MatrixDouble> jac,
4127 boost::shared_ptr<Range> edges_ptr)
4136 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
4139 return OP::doWork(side,
type, data);
4144 boost::shared_ptr<Range> edgesPtr;
4147 if (!geom_field_name.empty()) {
4148 auto jac = boost::make_shared<MatrixDouble>();
4149 auto det = boost::make_shared<VectorDouble>();
4151 geom_field_name, jac,
4153 : boost::make_shared<
Range>()));
4228 const auto nb_gauss_pts = getGaussPts().size2();
4230 dataAtPts->faceMaterialForceAtPts, nb_gauss_pts);
4231 dataAtPts->normalPressureAtPts.resize(nb_gauss_pts,
false);
4232 if (getNinTheLoop() == 0) {
4233 dataAtPts->faceMaterialForceAtPts.clear();
4236 auto loop_size = getLoopSize();
4237 if (loop_size == 1) {
4238 auto numebered_fe_ptr = getSidePtrFE()->numeredEntFiniteElementPtr;
4239 auto pstatus = numebered_fe_ptr->getPStatus();
4240 if (pstatus & (PSTATUS_SHARED | PSTATUS_MULTISHARED)) {
4247 auto t_normal = getFTensor1NormalsAtGaussPts();
4248 auto t_T =
dataAtPts->getFTensorFaceMaterialForce(
4252 auto t_P =
dataAtPts->getFTensorApproxP(nb_gauss_pts);
4253 auto t_u_gamma =
dataAtPts->getFTensorSmallHybridDisp(nb_gauss_pts);
4254 auto t_grad_u_gamma =
dataAtPts->getFTensorGradHybridDisp(nb_gauss_pts);
4255 auto t_strain =
dataAtPts->getFTensorLogStretch(nb_gauss_pts);
4256 auto t_omega =
dataAtPts->getFTensorRotAxis(nb_gauss_pts);
4277 case GRIFFITH_FORCE:
4278 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
4279 t_N(
I) = t_normal(
I);
4282 t_A(
i,
j) = levi_civita(
i,
j,
k) * t_omega(
k);
4284 t_grad_u(
i,
j) = t_R(
i,
j) + t_strain(
i,
j);
4286 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
4289 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
4292 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
4298 case GRIFFITH_SKELETON:
4299 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
4302 t_N(
I) = t_normal(
I);
4307 t_strain(
i,
j) - 0.5 * (t_grad_u_gamma(
i,
j) + t_grad_u_gamma(
j,
i));
4311 t_grad_u_gamma(
i,
J) +
4312 (2 * t_R(
i,
K) * t_N(
K) - (t_R(
k,
L) * t_N(
k) * t_N(
L)) * t_N(
i)) *
4315 t_T(
I) += t_N(
J) * (t_grad_u(
i,
I) * t_P(
i,
J)) / loop_size;
4318 t_T(
I) -= t_N(
I) * ((t_strain(
i,
K) * t_P(
i,
K)) / 2.) / loop_size;
4322 (t_N(
J) * ((
t_kd(
i,
I) + t_grad_u_gamma(
i,
I)) * t_P(
i,
J))) /
4331 "Grffith energy release "
4332 "selector not implemented");
4336 auto side_fe_ptr = getSidePtrFE();
4337 auto side_fe_mi_ptr = side_fe_ptr->numeredEntFiniteElementPtr;
4338 auto pstatus = side_fe_mi_ptr->getPStatus();
4340 auto owner = side_fe_mi_ptr->getOwnerProc();
4342 <<
"OpFaceSideMaterialForce: owner proc is not 0, owner proc: " << owner
4343 <<
" " << getPtrFE()->mField.get_comm_rank() <<
" n in the loop "
4344 << getNinTheLoop() <<
" loop size " << getLoopSize();
4356 auto fe_mi_ptr = getFEMethod()->numeredEntFiniteElementPtr;
4357 auto pstatus = fe_mi_ptr->getPStatus();
4359 auto owner = fe_mi_ptr->getOwnerProc();
4361 <<
"OpFaceMaterialForce: owner proc is not 0, owner proc: " << owner
4362 <<
" " << getPtrFE()->mField.get_comm_rank();
4370 double face_pressure = 0.;
4371 auto t_T =
dataAtPts->getFTensorFaceMaterialForce(
4372 getGaussPts().size2());
4375 auto t_w = getFTensor0IntegrationWeight();
4376 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
4377 t_face_T(
I) += t_w * t_T(
I);
4378 face_pressure += t_w * t_p;
4383 t_face_T(
I) *= getMeasure();
4384 face_pressure *= getMeasure();
4386 auto get_tag = [&](
auto name,
auto dim) {
4387 auto &moab = getPtrFE()->mField.get_moab();
4389 double def_val[] = {0., 0., 0.};
4390 CHK_MOAB_THROW(moab.tag_get_handle(name, dim, MB_TYPE_DOUBLE, tag,
4391 MB_TAG_CREAT | MB_TAG_SPARSE, def_val),
4396 auto set_tag = [&](
auto &&tag,
auto ptr) {
4397 auto &moab = getPtrFE()->mField.get_moab();
4398 auto face = getPtrFE()->getFEEntityHandle();
4399 CHK_MOAB_THROW(moab.tag_set_data(tag, &face, 1, ptr),
"set tag");
4402 set_tag(get_tag(
"MaterialForce", 3), &t_face_T(0));
4403 set_tag(get_tag(
"FacePressure", 1), &face_pressure);
4408template <
typename OP_PTR>
4409std::tuple<std::string, MatrixDouble>
4411 const std::string block_name) {
4413 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
4415 auto ts_time = op_ptr->getTStime();
4416 auto ts_time_step = op_ptr->getTStimeStep();
4423 MatrixDouble m_ref_coords = op_ptr->getCoordsAtGaussPts();
4424 MatrixDouble m_ref_normals = op_ptr->getNormalsAtGaussPts();
4426 auto v_analytical_expr =
4428 m_ref_coords, m_ref_normals, block_name);
4430 if (PetscUnlikely(!v_analytical_expr.size2())) {
4432 "Analytical expression is empty or does not exist, "
4433 "check python file");
4436 return std::make_tuple(block_name, v_analytical_expr);
4440 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4441 boost::shared_ptr<MatrixDouble> vec,
ScalarFun beta_coeff,
4442 boost::shared_ptr<Range> ents_ptr)
4443 :
OP(broken_base_side_data, ents_ptr) {
4444 this->sourceVec = vec;
4445 this->betaCoeff = beta_coeff;
4449 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4450 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
4451 :
OP(broken_base_side_data, ents_ptr) {
4452 this->sourceVec = boost::shared_ptr<MatrixDouble>();
4453 this->betaCoeff = beta_coeff;
4462 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
4467 if (!brokenBaseSideData) {
4472 auto do_work_rhs = [
this](
int row_side, EntityType row_type,
4480 OP::nbIntegrationPts = OP::getGaussPts().size2();
4482 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
4484 OP::locF.resize(OP::nbRows,
false);
4487 CHKERR this->iNtegrate(row_data);
4489 CHKERR this->aSsemble(row_data);
4493 switch (OP::opType) {
4495 for (
auto &bd : *brokenBaseSideData) {
4497 boost::shared_ptr<MatrixDouble>(brokenBaseSideData, &bd.getFlux());
4498 CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
4499 this->sourceVec.reset();
4504 (std::string(
"wrong op type ") +
4505 OpBaseDerivativesBase::OpTypeNames[OP::opType])
4513 const std::string row_field,
4514 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4515 ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr)
4516 :
OP(row_field, boost::shared_ptr<
MatrixDouble>(), beta_coeff, ents_ptr),
4517 brokenBaseSideDataPtr(broken_base_side_data) {
4518 this->betaCoeff = beta_coeff;
4524 for (
auto &bd : (*brokenBaseSideDataPtr)) {
4529 if (this->sourceVec->size2() !=
SPACE_DIM) {
4531 "Inconsistent size of the source vector");
4533 if (this->sourceVec->size1() != OP::getGaussPts().size2()) {
4535 "Inconsistent size of the source vector");
4539 CHKERR OP::iNtegrate(data);
4541 this->sourceVec.reset();
4547 std::string row_field,
4548 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4549 ScalarFun beta,
const bool assmb_transpose,
const bool only_transpose,
4550 boost::shared_ptr<Range> ents_ptr)
4551 :
OP(row_field, broken_base_side_data, assmb_transpose, only_transpose,
4553 this->betaCoeff = beta;
4558 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
4559 ScalarFun beta, boost::shared_ptr<Range> ents_ptr)
4560 :
OP(broken_base_side_data, ents_ptr) {
4562 this->betaCoeff = beta;
4563 OP::assembleTranspose =
false;
4564 OP::onlyTranspose =
false;
4573 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
4578 if (!brokenBaseSideData) {
4583 auto do_work_lhs = [
this](
int row_side,
int col_side, EntityType row_type,
4584 EntityType col_type,
4589 auto check_if_assemble_transpose = [&] {
4591 if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
4595 }
else if (OP::assembleTranspose) {
4601 OP::rowSide = row_side;
4602 OP::rowType = row_type;
4603 OP::colSide = col_side;
4604 OP::colType = col_type;
4606 OP::locMat.resize(OP::nbRows, OP::nbCols,
false);
4608 CHKERR this->iNtegrate(row_data, col_data);
4609 CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
4613 switch (OP::opType) {
4616 for (
auto &bd : *brokenBaseSideData) {
4619 if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
4621 "base functions not set");
4625 OP::nbRows = bd.getData().getIndices().size();
4628 OP::nbIntegrationPts = OP::getGaussPts().size2();
4629 OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(bd.getData());
4637 bd.getSide(), bd.getSide(),
4640 bd.getType(), bd.getType(),
4643 bd.getData(), bd.getData()
4652 (std::string(
"wrong op type ") +
4653 OpBaseDerivativesBase::OpTypeNames[OP::opType])
Auxilary functions for Eshelbian plasticity.
Eshelbian plasticity interface.
Lie algebra implementation.
#define FTENSOR_INDEXES(DIM,...)
#define FTENSOR_INDEX(DIM, I)
Fourth-order symmetrization tensor.
Fourth-order differential tensor symmetric in both index pairs.
Kronecker Delta class symmetric.
Mapping from symmetric tensor indices to packed storage index.
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 getDiffSpatialGradientDP(TInvD &t_d_u_d_b, TRotation &t_R)
std::tuple< std::string, MatrixDouble > getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr, const std::string block_name)
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)
auto getDiffSpatialGradientDR(TInvD &t_d_u_d_b, TRotation &t_R, TDiffRotation &t_diff_R, TStress &t_P)
boost::shared_ptr< VectorDouble > VectorPtr
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
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 StretchHandling stretchHandling
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static PetscBool physicalTimeFlg
static double currentPhysicalTime
static constexpr enum SymmetrySelector symmetrySelector
static boost::function< double(const double)> f
static PetscBool setSingularity
static bool hasNonHomogeneousMaterialBlock
static boost::function< double(const double)> d_f
static bool isNoStretch()
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.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
TSMethod::TSContext getTSCtx() const
double getTStimeStep() const
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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.
static constexpr Switches CtxSetTime
Time value switch.
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
double getVolume() const
element volume (linear geometry)
MoFEMErrorCode iNtegrate(EntData &data)
MatrixDouble K
local tangent matrix
VectorDouble nF
local right hand side vector
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode assemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
boost::shared_ptr< AnalyticalTractionBcVec > bcData
boost::shared_ptr< double > piolaScalePtr
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 &row_data, EntData &col_data)
boost::shared_ptr< double > piolaScalePtr
MoFEMErrorCode iNtegrate(EntData &data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< PressureBcVec > bcData
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
boost::shared_ptr< double > piolaScalePtr
boost::shared_ptr< TractionBcVec > bcData
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.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< ExternalStrainVec > externalStrainVecPtr
VectorPtr externalPressurePtr
OpCalculateExternalPressure(VectorPtr external_pressure_ptr, boost::shared_ptr< ExternalStrainVec > external_strain_vec_ptr, std::map< std::string, boost::shared_ptr< ScalingMethod > > smv)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
std::array< double, 6 > & reactionVec
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Caluclate face material force and normal pressure at gauss points.
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
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 &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &data)
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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)
const bool hasNonhomogeneousMatBlock
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
const bool hasNonhomogeneousMatBlock
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data)
const bool hasNonhomogeneousMatBlock
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)
double alphaViscousOmega0
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
double alphaViscousOmega0
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 &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
std::map< std::string, boost::shared_ptr< ScalingMethod > > scalingMethodsMap
MoFEMErrorCode iNtegrate(EntData &data)
