13 #include <boost/math/constants/constants.hpp>
38 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
41 t_R(
i,
j) += t_Omega(
i,
j);
46 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
47 const auto a = sin(angle) / angle;
48 t_R(
i,
j) +=
a * t_Omega(
i,
j);
52 const auto ss_2 = sin(angle / 2.) / angle;
53 const auto b = 2. * ss_2 * ss_2;
54 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
73 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
75 t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
79 const auto ss = sin(angle);
80 const auto a = ss / angle;
81 t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
84 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
86 const auto angle2 = angle * angle;
87 const auto cc = cos(angle);
88 const auto diff_a = (angle * cc - ss) / angle2;
89 t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
93 const auto ss_2 = sin(angle / 2.);
94 const auto cc_2 = cos(angle / 2.);
95 const auto b = 2. * ss_2 * ss_2 / angle2;
97 b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
98 FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
100 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
102 diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
107 template <
typename T>
116 constexpr
double eps = 1e-10;
117 for (
int l = 0;
l != 3; ++
l) {
119 t_omega_c(
i) = t_omega(
i);
120 t_omega_c(
l) += std::complex<double>(0,
eps);
122 for (
int i = 0;
i != 3; ++
i) {
123 for (
int j = 0;
j != 3; ++
j) {
124 for (
int k = 0;
k != 3; ++
k) {
125 t_diff2_R(
i,
j,
k,
l) = t_diff_R_c(
i,
j,
k).imag() /
eps;
135 static inline auto check(
const double &
a,
const double &b) {
136 constexpr
double eps = std::numeric_limits<float>::epsilon();
137 return std::abs(
a - b) / absMax <
eps;
142 double isEq::absMax = 1;
144 inline auto is_eq(
const double &
a,
const double &b) {
145 return isEq::check(
a, b);
149 std::array<double, DIM> tmp;
150 std::copy(ptr, &ptr[DIM], tmp.begin());
151 std::sort(tmp.begin(), tmp.end());
152 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
153 constexpr
double eps = std::numeric_limits<float>::epsilon();
154 isEq::absMax = std::max(isEq::absMax,
static_cast<double>(
eps));
155 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(),
is_eq));
163 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
165 int i = 0,
j = 1,
k = 2;
167 if (
is_eq(eig(0), eig(1))) {
171 }
else if (
is_eq(eig(0), eig(2))) {
175 }
else if (
is_eq(eig(1), eig(2))) {
182 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
184 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
186 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
194 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
207 int nb_integration_pts = getGaussPts().size2();
209 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
210 auto t_F = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
215 auto t_eshelby_stress =
216 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->SigmaAtPts);
220 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
221 t_eshelby_stress(
i,
j) = t_energy *
t_kd(
i,
j) - t_F(
m,
i) * t_P(
m,
j);
239 int nb_integration_pts = getGaussPts().size2();
247 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
248 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
249 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
250 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
251 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
252 dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts,
false);
254 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
255 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
256 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
257 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
259 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
260 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
261 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
262 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
263 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
264 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
265 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
267 dataAtPts->logStretch2H1AtPts.resize(6, nb_integration_pts,
false);
268 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
270 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
271 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
273 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
274 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
275 auto t_levi_kirchoff_dP =
276 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
277 auto t_levi_kirchoff_domega =
278 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
279 auto t_approx_P_adjont_dstretch =
280 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
281 auto t_approx_P_adjont_log_du =
282 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
283 auto t_approx_P_adjont_log_du_dP =
284 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
285 auto t_approx_P_adjont_log_du_domega =
286 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
287 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
288 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
289 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
291 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
292 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
293 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
295 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
296 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
297 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
299 auto &nbUniq = dataAtPts->nbUniq;
300 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
301 auto t_log_stretch_total =
302 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
304 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretch2H1AtPts);
307 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
316 switch (EshelbianCore::gradApproximator) {
329 auto calculate_rotation = [&]() {
333 t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
334 t_R(
i,
j) = t0(
i,
j);
337 auto calculate_stretch = [&]() {
340 t_u(
i,
j) = t_log_u(
i,
j) + t_kd_sym(
i,
j);
341 t_diff_u(
i,
j,
k,
l) = (t_kd_sym(
i,
k) ^ t_kd_sym(
j,
l)) / 4.;
342 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
344 eigen_vec(
i,
j) = t_log_u(
i,
j);
347 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
348 if (nbUniq[gg] < 3) {
349 sort_eigen_vals<3>(eig, eigen_vec);
351 t_eigen_vals(
i) = eig(
i);
352 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
355 t_u(
i,
j) = t_u_tmp(
i,
j);
359 t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
360 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
364 calculate_rotation();
366 if (!EshelbianCore::noStretch) {
370 switch (EshelbianCore::gradApproximator) {
375 t_Ru(
i,
m) = t_R(
i,
l) * t_u(
l,
m);
376 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
377 t_h_domega(
i,
j,
k) = (t_diff_R(
i,
l,
k) * t_u(
l,
m)) * t_h1(
m,
j);
378 t_h_dlog_u(
i,
j,
L) = (t_R(
i,
l) * t_Ldiff_u(
l,
m,
L)) * t_h1(
m,
j);
381 t_approx_P_intermidiate(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
382 t_approx_P_adjont_dstretch(
l,
m) =
383 t_approx_P_intermidiate(
i,
m) * t_R(
i,
l);
387 t_levi_kirchoff_dP(
i,
j,
k) =
389 t_levi_kirchoff_domega(
k,
n) =
391 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
393 t_approx_P_adjont_log_du(
L) =
394 t_Ldiff_u(
l,
m,
L) * t_approx_P_adjont_dstretch(
l,
m);
395 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
396 t_approx_P_adjont_log_du_domega(
n,
L) =
398 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
407 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
408 t_Ru(
i,
m) = t_Omega(
i,
m) + t_u(
i,
m);
409 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
410 t_h_domega(
i,
j,
k) =
411 FTensor::levi_civita<double>(
i,
m,
k) * t_h1(
m,
j);
412 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
m,
L) * t_h1(
m,
j);
415 t_approx_P_intermidiate(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
416 t_approx_P_adjont_dstretch(
i,
m) = t_approx_P_intermidiate(
i,
m);
421 t_levi_kirchoff_domega(
k,
n) = 0;
423 t_approx_P_adjont_log_du(
L) =
424 t_Ldiff_u(
i,
m,
L) * t_approx_P_adjont_dstretch(
i,
m);
425 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
426 t_approx_P_adjont_log_du_domega(
n,
L) = 0;
438 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
439 t_h(
i,
j) = t_Omega(
i,
j) + t_u(
i,
j);
440 t_h_domega(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
441 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
446 t_levi_kirchoff_domega(
k,
l) = 0;
448 t_approx_P_adjont_dstretch(
i,
j) = t_approx_P(
i,
j);
449 t_approx_P_adjont_log_du(
L) =
450 t_Ldiff_u(
i,
j,
L) * t_approx_P_adjont_dstretch(
i,
j);
451 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
452 t_approx_P_adjont_log_du_domega(
m,
L) = 0;
460 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
461 eigen_vec(
i,
j) = t_C_h1(
i,
j);
463 switch (EshelbianCore::stretchSelector) {
475 switch (EshelbianCore::gradApproximator) {
477 t_log_u2_h1(
i,
j) = 0;
478 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
482 auto t_log_u2_h1_tmp =
484 t_log_u2_h1(
i,
j) = t_log_u2_h1_tmp(
i,
j);
485 t_log_stretch_total(
i,
j) = t_log_u2_h1_tmp(
i,
j) / 2 + t_log_u(
i,
j);
488 t_log_u2_h1(
i,
j) = 0;
489 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
498 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
500 t_h(
i,
j) = t_Omega(
i,
j) + t_u(
i,
j);
501 t_h_domega(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
506 t_levi_kirchoff_domega(
k,
l) = 0;
508 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
515 ++t_levi_kirchoff_dP;
516 ++t_levi_kirchoff_domega;
517 ++t_approx_P_adjont_dstretch;
518 ++t_approx_P_adjont_log_du;
519 ++t_approx_P_adjont_log_du_dP;
520 ++t_approx_P_adjont_log_du_domega;
532 ++t_log_stretch_total;
541 int nb_integration_pts = data.
getN().size1();
542 auto v = getVolume();
543 auto t_w = getFTensor0IntegrationWeight();
544 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
545 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
546 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
547 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
548 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
549 dataAtPts->wL2DotDotAtPts.clear();
551 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
553 int nb_base_functions = data.
getN().size2();
557 auto get_ftensor1 = [](
auto &
v) {
562 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
564 auto t_nf = get_ftensor1(nF);
566 for (; bb != nb_dofs / 3; ++bb) {
567 t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
568 t_nf(
i) -=
a * t_row_base_fun * alphaW * t_s_dot_w(
i);
569 t_nf(
i) -=
a * t_row_base_fun * alphaRho * t_s_dot_dot_w(
i);
573 for (; bb != nb_base_functions; ++bb)
587 int nb_integration_pts = getGaussPts().size2();
588 auto v = getVolume();
589 auto t_w = getFTensor0IntegrationWeight();
590 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
591 int nb_base_functions = data.
getN().size2();
596 auto get_ftensor1 = [](
auto &
v) {
600 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
602 auto t_nf = get_ftensor1(nF);
604 for (; bb != nb_dofs / 3; ++bb) {
605 t_nf(
k) += (
a * t_row_base_fun) * t_levi_kirchoff(
k);
609 for (; bb != nb_base_functions; ++bb) {
621 int nb_integration_pts = data.
getN().size1();
622 auto v = getVolume();
623 auto t_w = getFTensor0IntegrationWeight();
624 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
626 int nb_base_functions = data.
getN().size2() / 3;
633 auto get_ftensor1 = [](
auto &
v) {
638 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
640 auto t_nf = get_ftensor1(nF);
647 for (; bb != nb_dofs / 3; ++bb) {
648 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
653 for (; bb != nb_base_functions; ++bb)
666 int nb_integration_pts = data.
getN().size1();
667 auto v = getVolume();
668 auto t_w = getFTensor0IntegrationWeight();
669 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
671 int nb_base_functions = data.
getN().size2() / 9;
678 auto get_ftensor0 = [](
auto &
v) {
682 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
684 auto t_nf = get_ftensor0(nF);
691 for (; bb != nb_dofs; ++bb) {
692 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
696 for (; bb != nb_base_functions; ++bb) {
709 int nb_integration_pts = data.
getN().size1();
710 auto v = getVolume();
711 auto t_w = getFTensor0IntegrationWeight();
712 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
713 int nb_base_functions = data.
getN().size2() / 3;
716 auto get_ftensor1 = [](
auto &
v) {
721 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
723 auto t_nf = get_ftensor1(nF);
725 for (; bb != nb_dofs / 3; ++bb) {
726 double div_row_base = t_row_diff_base_fun(
i,
i);
727 t_nf(
i) +=
a * div_row_base * t_w_l2(
i);
729 ++t_row_diff_base_fun;
731 for (; bb != nb_base_functions; ++bb) {
732 ++t_row_diff_base_fun;
741 template <AssemblyType A>
747 for (
auto &bc : (*bcDispPtr)) {
749 if (bc.faces.find(fe_ent) != bc.faces.end()) {
751 int nb_integration_pts = OP::getGaussPts().size2();
752 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
753 auto t_w = OP::getFTensor0IntegrationWeight();
754 int nb_base_functions = data.
getN().size2() / 3;
761 for (
auto &sm : scalingMethodsVec) {
762 scale *= sm->getScale(OP::getFEMethod()->ts_t);
769 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
770 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
772 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
774 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_disp(
i) * 0.5;
778 for (; bb != nb_base_functions; ++bb)
790 return OP::iNtegrate(data);
793 template <AssemblyType A>
804 for (
auto &bc : (*bcRotPtr)) {
806 if (bc.faces.find(fe_ent) != bc.faces.end()) {
808 int nb_integration_pts = OP::getGaussPts().size2();
809 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
810 auto t_w = OP::getFTensor0IntegrationWeight();
812 int nb_base_functions = data.
getN().size2() / 3;
815 auto get_ftensor1 = [](
auto &
v) {
823 double theta = bc.theta;
824 theta *= OP::getFEMethod()->ts_t;
827 const double a = sqrt(t_normal(
i) * t_normal(
i));
828 t_omega(
i) = t_normal(
i) * (theta /
a);
830 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
832 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
834 t_delta(
i) = t_center(
i) - t_coords(
i);
836 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
838 auto t_nf = getFTensor1FromPtr<3>(&*OP::locF.begin());
840 for (; bb != nb_dofs / 3; ++bb) {
841 t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
845 for (; bb != nb_base_functions; ++bb)
858 return OP::iNtegrate(data);
867 int nb_integration_pts = getGaussPts().size2();
868 int nb_base_functions = data.
getN().size2();
869 double ts_t = getFEMethod()->ts_t;
872 if (this->locF.size() != nb_dofs)
874 "Size of locF %d != nb_dofs %d", this->locF.size(), nb_dofs);
877 auto integrate_rhs = [&](
auto &bc) {
880 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
882 auto t_w = getFTensor0IntegrationWeight();
884 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
885 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
887 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
888 t_f(
i) -= ts_t * t_w * t_row_base * t_val(
i);
892 for (; rr != nb_base_functions; ++rr)
897 this->locF *= getMeasure();
901 auto integrate_rhs_cook = [&](
auto &bc) {
904 auto t_val = getFTensor1FromPtr<3>(&*bc.vals.begin());
906 auto t_w = getFTensor0IntegrationWeight();
907 auto t_coords = getFTensor1CoordsAtGaussPts();
909 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
911 auto calc_tau = [](
double y) {
914 return -y * (y - 1) / 0.25;
917 const auto tau = calc_tau(t_coords(1));
918 auto t_f = getFTensor1FromPtr<3>(&*this->locF.begin());
920 for (; rr != nb_dofs /
SPACE_DIM; ++rr) {
921 t_f(
i) -= ts_t * t_w * t_row_base * tau * t_val(
i);
926 for (; rr != nb_base_functions; ++rr)
934 this->locF *= 2. * getMeasure();
941 for (
auto &bc : *(bcData)) {
942 if (bc.faces.find(fe_ent) != bc.faces.end()) {
947 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*")))
948 CHKERR integrate_rhs_cook(bc);
961 int nb_integration_pts = row_data.
getN().size1();
962 int row_nb_dofs = row_data.
getIndices().size();
963 int col_nb_dofs = col_data.
getIndices().size();
966 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2));
969 auto v = getVolume();
970 auto t_w = getFTensor0IntegrationWeight();
971 int row_nb_base_functions = row_data.
getN().size2();
973 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
976 for (; rr != row_nb_dofs / 3; ++rr) {
978 auto t_m = get_ftensor1(
K, 3 * rr, 0);
979 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
980 double div_col_base = t_col_diff_base_fun(
i,
i);
981 t_m(
i) +=
a * t_row_base_fun * div_col_base;
983 ++t_col_diff_base_fun;
987 for (; rr != row_nb_base_functions; ++rr)
998 if (alphaW < std::numeric_limits<double>::epsilon() &&
999 alphaRho < std::numeric_limits<double>::epsilon())
1002 const int nb_integration_pts = row_data.
getN().size1();
1003 const int row_nb_dofs = row_data.
getIndices().size();
1006 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2)
1012 auto v = getVolume();
1013 auto t_w = getFTensor0IntegrationWeight();
1015 int row_nb_base_functions = row_data.
getN().size2();
1018 double ts_scale = alphaW * getTSa();
1019 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1020 ts_scale += alphaRho * getTSaa();
1022 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1023 double a =
v * t_w * ts_scale;
1026 for (; rr != row_nb_dofs / 3; ++rr) {
1029 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1030 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1031 const double b =
a * t_row_base_fun * t_col_base_fun;
1040 for (; rr != row_nb_base_functions; ++rr)
1055 int nb_integration_pts = row_data.
getN().size1();
1056 int row_nb_dofs = row_data.
getIndices().size();
1057 int col_nb_dofs = col_data.
getIndices().size();
1061 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1063 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1065 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
1067 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
1069 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
1071 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2));
1077 auto v = getVolume();
1078 auto t_w = getFTensor0IntegrationWeight();
1079 auto t_approx_P_adjont_log_du_dP =
1080 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1081 int row_nb_base_functions = row_data.
getN().size2();
1083 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1086 for (; rr != row_nb_dofs / 6; ++rr) {
1089 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1091 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1093 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
1101 for (; rr != row_nb_base_functions; ++rr)
1104 ++t_approx_P_adjont_log_du_dP;
1115 int nb_integration_pts = row_data.
getN().size1();
1116 int row_nb_dofs = row_data.
getIndices().size();
1117 int col_nb_dofs = col_data.
getIndices().size();
1120 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c), &
m(
r + 3,
c), &
m(
r + 4,
c),
1127 auto v = getVolume();
1128 auto t_w = getFTensor0IntegrationWeight();
1129 auto t_approx_P_adjont_log_du_dP =
1130 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1132 int row_nb_base_functions = row_data.
getN().size2();
1134 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1137 for (; rr != row_nb_dofs / 6; ++rr) {
1138 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1139 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1140 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1142 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
1149 for (; rr != row_nb_base_functions; ++rr)
1152 ++t_approx_P_adjont_log_du_dP;
1164 int row_nb_dofs = row_data.
getIndices().size();
1165 int col_nb_dofs = col_data.
getIndices().size();
1169 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1171 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1173 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
1175 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
1177 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
1179 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2)
1189 auto v = getVolume();
1190 auto t_w = getFTensor0IntegrationWeight();
1191 auto t_approx_P_adjont_log_du_domega =
1192 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
1194 int row_nb_base_functions = row_data.
getN().size2();
1197 int nb_integration_pts = row_data.
getN().size1();
1199 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1203 for (; rr != row_nb_dofs / 6; ++rr) {
1205 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1206 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1207 double v =
a * t_row_base_fun * t_col_base_fun;
1208 t_m(
L,
k) +=
v * t_approx_P_adjont_log_du_domega(
k,
L);
1215 for (; rr != row_nb_base_functions; ++rr)
1219 ++t_approx_P_adjont_log_du_domega;
1228 int nb_integration_pts = getGaussPts().size2();
1229 int row_nb_dofs = row_data.
getIndices().size();
1230 int col_nb_dofs = col_data.
getIndices().size();
1234 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1236 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1238 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1248 auto v = getVolume();
1249 auto t_w = getFTensor0IntegrationWeight();
1250 auto t_levi_kirchoff_dP =
1251 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1253 int row_nb_base_functions = row_data.
getN().size2();
1256 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1259 for (; rr != row_nb_dofs / 3; ++rr) {
1260 double b =
a * t_row_base_fun;
1262 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1263 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1264 t_m(
k,
i) += b * (t_levi_kirchoff_dP(
i,
l,
k) * t_col_base_fun(
l));
1270 for (; rr != row_nb_base_functions; ++rr) {
1275 ++t_levi_kirchoff_dP;
1283 int nb_integration_pts = getGaussPts().size2();
1284 int row_nb_dofs = row_data.
getIndices().size();
1285 int col_nb_dofs = col_data.
getIndices().size();
1289 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c));
1298 auto v = getVolume();
1299 auto t_w = getFTensor0IntegrationWeight();
1300 auto t_levi_kirchoff_dP =
1301 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
1303 int row_nb_base_functions = row_data.
getN().size2();
1306 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1309 for (; rr != row_nb_dofs / 3; ++rr) {
1310 double b =
a * t_row_base_fun;
1311 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1312 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1313 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1314 t_m(
k) += b * (t_levi_kirchoff_dP(
i,
j,
k) * t_col_base_fun(
i,
j));
1321 for (; rr != row_nb_base_functions; ++rr) {
1325 ++t_levi_kirchoff_dP;
1333 int nb_integration_pts = getGaussPts().size2();
1334 int row_nb_dofs = row_data.
getIndices().size();
1335 int col_nb_dofs = col_data.
getIndices().size();
1339 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1341 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1343 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1354 auto v = getVolume();
1355 auto t_w = getFTensor0IntegrationWeight();
1356 auto t_levi_kirchoff_domega =
1357 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
1358 int row_nb_base_functions = row_data.
getN().size2();
1360 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1363 for (; rr != row_nb_dofs / 3; ++rr) {
1364 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1365 const double b =
a * t_row_base_fun;
1367 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1368 t_m(
k,
l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(
k,
l);
1374 for (; rr != row_nb_base_functions; ++rr) {
1378 ++t_levi_kirchoff_domega;
1386 if (dataAtPts->matInvD.size1() ==
size_symm &&
1387 dataAtPts->matInvD.size2() ==
size_symm) {
1388 return integrateImpl<0>(row_data, col_data);
1390 return integrateImpl<size_symm>(row_data, col_data);
1403 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1405 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1407 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1412 int nb_integration_pts = getGaussPts().size2();
1413 int row_nb_dofs = row_data.
getIndices().size();
1414 int col_nb_dofs = col_data.
getIndices().size();
1416 auto v = getVolume();
1417 auto t_w = getFTensor0IntegrationWeight();
1418 int row_nb_base_functions = row_data.
getN().size2() / 3;
1425 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1426 &*dataAtPts->matInvD.data().begin());
1429 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1433 for (; rr != row_nb_dofs / 3; ++rr) {
1435 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1436 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1437 t_m(
i,
k) +=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
1445 for (; rr != row_nb_base_functions; ++rr)
1455 OpSpatialConsistency_dBubble_dBubble::integrate(
EntData &row_data,
1458 if (dataAtPts->matInvD.size1() ==
size_symm &&
1459 dataAtPts->matInvD.size2() ==
size_symm) {
1460 return integrateImpl<0>(row_data, col_data);
1462 return integrateImpl<size_symm>(row_data, col_data);
1469 OpSpatialConsistency_dBubble_dBubble::integrateImpl(
EntData &row_data,
1473 int nb_integration_pts = getGaussPts().size2();
1474 int row_nb_dofs = row_data.
getIndices().size();
1475 int col_nb_dofs = col_data.
getIndices().size();
1477 auto v = getVolume();
1478 auto t_w = getFTensor0IntegrationWeight();
1479 int row_nb_base_functions = row_data.
getN().size2() / 9;
1486 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1487 &*dataAtPts->matInvD.data().begin());
1490 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1494 for (; rr != row_nb_dofs; ++rr) {
1496 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1498 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
1505 for (; rr != row_nb_base_functions; ++rr)
1516 if (dataAtPts->matInvD.size1() ==
size_symm &&
1517 dataAtPts->matInvD.size2() ==
size_symm) {
1518 return integrateImpl<0>(row_data, col_data);
1520 return integrateImpl<size_symm>(row_data, col_data);
1533 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2)
1538 int nb_integration_pts = getGaussPts().size2();
1539 int row_nb_dofs = row_data.
getIndices().size();
1540 int col_nb_dofs = col_data.
getIndices().size();
1542 auto v = getVolume();
1543 auto t_w = getFTensor0IntegrationWeight();
1544 int row_nb_base_functions = row_data.
getN().size2() / 9;
1553 auto t_inv_D = getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, S>(
1554 &*dataAtPts->matInvD.data().begin());
1557 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1560 auto t_m = get_ftensor1(
K, 0, 0);
1563 for (; rr != row_nb_dofs; ++rr) {
1565 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1566 t_m(
k) +=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
1574 for (; rr != row_nb_base_functions; ++rr)
1586 int nb_integration_pts = row_data.
getN().size1();
1587 int row_nb_dofs = row_data.
getIndices().size();
1588 int col_nb_dofs = col_data.
getIndices().size();
1593 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1595 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1597 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
1609 auto v = getVolume();
1610 auto t_w = getFTensor0IntegrationWeight();
1611 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1612 int row_nb_base_functions = row_data.
getN().size2() / 3;
1615 const double ts_a = getTSa();
1617 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1623 for (; rr != row_nb_dofs / 3; ++rr) {
1626 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
1629 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1630 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1631 t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
1639 for (; rr != row_nb_base_functions; ++rr)
1648 OpSpatialConsistency_dBubble_domega::integrate(
EntData &row_data,
1652 int nb_integration_pts = row_data.
getN().size1();
1653 int row_nb_dofs = row_data.
getIndices().size();
1654 int col_nb_dofs = col_data.
getIndices().size();
1658 &
m(
r,
c + 0), &
m(
r,
c + 1), &
m(
r,
c + 2));
1668 auto v = getVolume();
1669 auto t_w = getFTensor0IntegrationWeight();
1670 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
1671 int row_nb_base_functions = row_data.
getN().size2() / 9;
1674 const double ts_a = getTSa();
1676 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1680 for (; rr != row_nb_dofs; ++rr) {
1683 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
1686 auto t_m = get_ftensor2(
K, rr, 0);
1687 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1688 t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
1696 for (; rr != row_nb_base_functions; ++rr)
1709 if(tagSense != getSkeletonSense())
1712 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1713 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1715 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
1716 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1721 Tag th_w = create_tag(
"SpatialDisplacement", 3);
1722 Tag th_omega = create_tag(
"Omega", 3);
1723 Tag th_approxP = create_tag(
"Piola1Stress", 9);
1724 Tag th_calcSigma = create_tag(
"EshelbyStress", 9);
1725 Tag th_sigma = create_tag(
"CauchyStress", 9);
1726 Tag th_log_u = create_tag(
"LogSpatialStretch", 9);
1727 Tag th_u = create_tag(
"SpatialStretch", 9);
1728 Tag th_h = create_tag(
"h", 9);
1729 Tag th_X = create_tag(
"X", 3);
1730 Tag th_detF = create_tag(
"detF", 1);
1731 Tag th_angular_momentum = create_tag(
"AngularMomentum", 3);
1735 Tag th_traction = create_tag(
"traction", 3);
1737 Tag th_disp = create_tag(
"H1Displacement", 3);
1738 Tag th_disp_error = create_tag(
"DisplacementError", 1);
1739 Tag th_lambda_disp = create_tag(
"ContactDisplacement", 3);
1741 Tag th_energy = create_tag(
"Energy", 1);
1743 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1744 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
1745 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1747 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
1748 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1749 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1750 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
1751 if (dataAtPts->SigmaAtPts.size2() != dataAtPts->approxPAtPts.size2()) {
1753 dataAtPts->SigmaAtPts.resize(dataAtPts->approxPAtPts.size1(),
1754 dataAtPts->approxPAtPts.size2(),
false);
1755 dataAtPts->SigmaAtPts.clear();
1758 auto t_calcSigma_P = getFTensor2FromMat<3, 3>(dataAtPts->SigmaAtPts);
1759 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
1760 auto t_coords = getFTensor1FromMat<3>(dataAtPts->XH1AtPts);
1761 auto t_normal = getFTensor1NormalsAtGaussPts();
1762 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
1763 auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
1765 if (dataAtPts->energyAtPts.size() == 0) {
1767 dataAtPts->energyAtPts.resize(getGaussPts().size2());
1768 dataAtPts->energyAtPts.clear();
1777 auto set_float_precision = [](
const double x) {
1778 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1785 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1787 v = set_float_precision(
v);
1788 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
1795 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1798 for (
auto &
a :
v.data())
1799 a = set_float_precision(
a);
1800 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
1801 &*
v.data().begin());
1809 &
m(0, 0), &
m(0, 1), &
m(0, 2),
1811 &
m(1, 0), &
m(1, 1), &
m(1, 2),
1813 &
m(2, 0), &
m(2, 1), &
m(2, 2));
1815 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1817 t_m(
i,
j) = t_d(
i,
j);
1818 for (
auto &
v :
m.data())
1819 v = set_float_precision(
v);
1820 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
1821 &*
m.data().begin());
1825 const auto nb_gauss_pts = getGaussPts().size2();
1826 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1829 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
1832 CHKERR save_vec_tag(th_w, t_w, gg);
1833 CHKERR save_vec_tag(th_X, t_coords, gg);
1834 CHKERR save_vec_tag(th_omega, t_omega, gg);
1835 CHKERR save_vec_tag(th_traction, t_traction, gg);
1838 CHKERR save_mat_tag(th_h, t_h, gg);
1841 for (
int ii = 0; ii != 3; ++ii)
1842 for (
int jj = 0; jj != 3; ++jj)
1843 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
1845 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
1848 for (
int ii = 0; ii != 3; ++ii)
1849 for (
int jj = 0; jj != 3; ++jj)
1850 t_u_tmp(ii, jj) = t_u(ii, jj);
1852 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
1853 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
1854 CHKERR save_mat_tag(th_calcSigma, t_calcSigma_P, gg);
1855 CHKERR save_vec_tag(th_disp, t_disp, gg);
1856 CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
1858 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
1859 CHKERR save_scal_tag(th_disp_error, u_error, gg);
1860 CHKERR save_scal_tag(th_energy, t_energy, gg);
1864 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
1865 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
1866 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
1869 t_levi(
k) = t_levi_kirchoff(
k);
1870 CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
1873 auto get_eiegn_vector_symmetric = [&](
auto &t_u) {
1893 CHKERR get_eiegn_vector_symmetric(t_u);
1915 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1916 std::vector<FieldSpace> spaces, std::string geom_field_name,
1917 boost::shared_ptr<Range> crack_front_edges_ptr) {
1928 boost::shared_ptr<Range> edges_ptr)
1937 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
1940 return OP::doWork(side,
type, data);
1945 boost::shared_ptr<Range> edgesPtr;
1948 auto jac = boost::make_shared<MatrixDouble>();
1949 auto det = boost::make_shared<VectorDouble>();
1951 geom_field_name, EshelbianCore::setSingularity
1952 ? crack_front_edges_ptr
1953 : boost::make_shared<Range>()));
1966 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
1967 std::vector<FieldSpace> spaces, std::string geom_field_name,
1968 boost::shared_ptr<Range> crack_front_edges_ptr) {
1972 auto jac = boost::make_shared<MatrixDouble>();
1973 auto det = boost::make_shared<VectorDouble>();
1975 geom_field_name, jac));
1989 boost::shared_ptr<MatrixDouble> jac,
1990 boost::shared_ptr<Range> edges_ptr)
1999 if (
type == MBEDGE && edgesPtr->find(ent) != edgesPtr->end()) {
2002 return OP::doWork(side,
type, data);
2007 boost::shared_ptr<Range> edgesPtr;
2010 auto jac = boost::make_shared<MatrixDouble>();
2011 auto det = boost::make_shared<VectorDouble>();
2013 geom_field_name, jac,
2014 EshelbianCore::setSingularity ? crack_front_edges_ptr
2015 : boost::make_shared<Range>()));
2022 nullptr,
nullptr,
nullptr);