12 #include <boost/math/constants/constants.hpp>
27 t_diff_deviator(
i,
j,
k,
l) = t_diff_stress(
i,
j,
k,
l);
29 constexpr
double third = boost::math::constants::third<double>();
31 t_diff_deviator(0, 0, 0, 0) -=
third;
32 t_diff_deviator(0, 0, 1, 1) -=
third;
34 t_diff_deviator(1, 1, 0, 0) -=
third;
35 t_diff_deviator(1, 1, 1, 1) -=
third;
37 t_diff_deviator(2, 2, 0, 0) -=
third;
38 t_diff_deviator(2, 2, 1, 1) -=
third;
40 t_diff_deviator(0, 0, 2, 2) -=
third;
41 t_diff_deviator(1, 1, 2, 2) -=
third;
42 t_diff_deviator(2, 2, 2, 2) -=
third;
44 return t_diff_deviator;
47 constexpr
static auto size_symm = (3 * (3 + 1)) / 2;
84 t_diff(
i,
j,
k,
l) = 0;
85 t_diff(0, 0, 0, 0) = 1;
86 t_diff(1, 1, 1, 1) = 1;
88 t_diff(1, 0, 1, 0) = 0.5;
89 t_diff(1, 0, 0, 1) = 0.5;
91 t_diff(0, 1, 0, 1) = 0.5;
92 t_diff(0, 1, 1, 0) = 0.5;
95 t_diff(2, 2, 2, 2) = 1;
97 t_diff(2, 0, 2, 0) = 0.5;
98 t_diff(2, 0, 0, 2) = 0.5;
99 t_diff(0, 2, 0, 2) = 0.5;
100 t_diff(0, 2, 2, 0) = 0.5;
102 t_diff(2, 1, 2, 1) = 0.5;
103 t_diff(2, 1, 1, 2) = 0.5;
104 t_diff(1, 2, 1, 2) = 0.5;
105 t_diff(1, 2, 2, 1) = 0.5;
118 template <
typename T>
133 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
136 t_R(
i,
j) += t_Omega(
i,
j);
141 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
142 const auto a = sin(angle) / angle;
143 t_R(
i,
j) +=
a * t_Omega(
i,
j);
147 const auto ss_2 = sin(angle / 2.) / angle;
148 const auto b = 2. * ss_2 * ss_2;
149 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
154 template <
typename T>
168 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
170 t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
174 const auto ss = sin(angle);
175 const auto a = ss / angle;
176 t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
179 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
181 const auto angle2 = angle * angle;
182 const auto cc = cos(angle);
183 const auto diff_a = (angle * cc - ss) / angle2;
184 t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
188 const auto ss_2 = sin(angle / 2.);
189 const auto cc_2 = cos(angle / 2.);
190 const auto b = 2. * ss_2 * ss_2 / angle2;
192 b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
193 FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
195 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
197 diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
202 template <
typename T>
211 constexpr
double eps = 1e-10;
212 for (
int l = 0;
l != 3; ++
l) {
214 t_omega_c(
i) = t_omega(
i);
215 t_omega_c(
l) += std::complex<double>(0,
eps);
217 for (
int i = 0;
i != 3; ++
i) {
218 for (
int j = 0;
j != 3; ++
j) {
219 for (
int k = 0;
k != 3; ++
k) {
220 t_diff2_R(
i,
j,
k,
l) = t_diff_R_c(
i,
j,
k).imag() /
eps;
230 static inline auto check(
const double &
a,
const double &b) {
231 constexpr
double eps = std::numeric_limits<float>::epsilon();
232 return std::abs(
a - b) / absMax <
eps;
237 double isEq::absMax = 1;
239 inline auto is_eq(
const double &
a,
const double &b) {
240 return isEq::check(
a, b);
244 std::array<double, DIM> tmp;
245 std::copy(ptr, &ptr[DIM], tmp.begin());
246 std::sort(tmp.begin(), tmp.end());
247 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
248 constexpr
double eps = std::numeric_limits<float>::epsilon();
249 isEq::absMax = std::max(isEq::absMax,
static_cast<double>(
eps));
250 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(),
is_eq));
258 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
260 int i = 0,
j = 1,
k = 2;
262 if (
is_eq(eig(0), eig(1))) {
266 }
else if (
is_eq(eig(0), eig(2))) {
270 }
else if (
is_eq(eig(1), eig(2))) {
277 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
279 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
281 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
289 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
304 auto nb_integration_pts = dataAtPts->approxPAtPts.size2();
306 if (nb_integration_pts != getGaussPts().size2()) {
308 "inconsistent number of integration points");
312 auto get_invert_D = [&]() {
333 dataAtPts->logStretchTensorAtPts.resize(
size_symm, nb_integration_pts,
false);
336 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
338 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
340 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD);
344 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
345 t_log_u(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_approx_P(
k,
l);
362 int nb_integration_pts = getGaussPts().size2();
370 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
371 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
372 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
373 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
374 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
375 dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts,
false);
377 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
378 dataAtPts->adjointPdUAtPts.resize(
size_symm, nb_integration_pts,
false);
379 dataAtPts->adjointPdUdPAtPts.resize(9 *
size_symm, nb_integration_pts,
false);
380 dataAtPts->adjointPdUdOmegaAtPts.resize(3 *
size_symm, nb_integration_pts,
382 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
383 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
384 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
385 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
386 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
387 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
388 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
390 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
392 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
393 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
395 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->hdLogStretchAtPts);
396 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
397 auto t_levi_kirchoff_dP =
398 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
399 auto t_levi_kirchoff_domega =
400 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
401 auto t_approx_P_adjont_dstretch =
402 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
403 auto t_approx_P_adjont_log_du =
404 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
405 auto t_approx_P_adjont_log_du_dP =
406 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
407 auto t_approx_P_adjont_log_du_domega =
408 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
409 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
410 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
411 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(dataAtPts->diffRotMatAtPts);
413 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
414 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
415 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
417 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
418 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
419 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
421 auto &nbUniq = dataAtPts->nbUniq;
422 auto t_grad_h1 = getFTensor2FromMat<3, 3>(dataAtPts->wGradH1AtPts);
423 auto t_log_stretch_total =
424 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
427 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
436 switch (EshelbianCore::gradApperoximator) {
449 auto calculate_rotation = [&]() {
453 t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
454 t_R(
i,
j) = t0(
i,
j);
457 auto calculate_stretch = [&]() {
458 eigen_vec(
i,
j) = t_log_u(
i,
j);
461 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
462 if (nbUniq[gg] < 3) {
463 sort_eigen_vals<3>(eig, eigen_vec);
465 t_eigen_vals(
i) = eig(
i);
466 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
469 t_u(
i,
j) = t_u_tmp(
i,
j);
473 t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
474 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
477 calculate_rotation();
479 if (!EshelbianCore::noStretch) {
483 switch (EshelbianCore::gradApperoximator) {
488 t_Ru(
i,
m) = t_R(
i,
l) * t_u(
l,
m);
489 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
490 t_h_domega(
i,
j,
k) = (t_diff_R(
i,
l,
k) * t_u(
l,
m)) * t_h1(
m,
j);
491 t_h_dlog_u(
i,
j,
L) = (t_R(
i,
l) * t_Ldiff_u(
l,
m,
L)) * t_h1(
m,
j);
494 t_approx_P_intermidiate(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
495 t_approx_P_adjont_dstretch(
l,
m) =
496 t_approx_P_intermidiate(
i,
m) * t_R(
i,
l);
500 t_levi_kirchoff_dP(
i,
j,
k) =
502 t_levi_kirchoff_domega(
k,
n) =
504 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
506 t_approx_P_adjont_log_du(
L) =
507 t_Ldiff_u(
l,
m,
L) * t_approx_P_adjont_dstretch(
l,
m);
508 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
509 t_approx_P_adjont_log_du_domega(
n,
L) =
511 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
520 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
521 t_Ru(
i,
m) = t_Omega(
i,
m) + t_u(
i,
m);
522 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
523 t_h_domega(
i,
j,
k) =
524 FTensor::levi_civita<double>(
i,
m,
k) * t_h1(
m,
j);
525 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
m,
L) * t_h1(
m,
j);
528 t_approx_P_intermidiate(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
529 t_approx_P_adjont_dstretch(
i,
m) = t_approx_P_intermidiate(
i,
m);
534 t_levi_kirchoff_domega(
k,
n) = 0;
536 t_approx_P_adjont_log_du(
L) =
537 t_Ldiff_u(
i,
m,
L) * t_approx_P_adjont_dstretch(
i,
m);
538 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
539 t_approx_P_adjont_log_du_domega(
n,
L) = 0;
551 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
552 t_h(
i,
j) = t_Omega(
i,
j) + t_u(
i,
j);
553 t_h_domega(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
554 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
559 t_levi_kirchoff_domega(
k,
l) = 0;
561 t_approx_P_adjont_dstretch(
i,
j) = t_approx_P(
i,
j);
562 t_approx_P_adjont_log_du(
L) =
563 t_Ldiff_u(
i,
j,
L) * t_approx_P_adjont_dstretch(
i,
j);
564 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
565 t_approx_P_adjont_log_du_domega(
m,
L) = 0;
573 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
574 eigen_vec(
i,
j) = t_C_h1(
i,
j);
576 switch (EshelbianCore::stretchSelector) {
580 for (
int ii = 0; ii != 3; ++ii) {
583 std::numeric_limits<float>::min_exponent)));
588 switch (EshelbianCore::gradApperoximator) {
590 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
594 auto t_log_u2_h1_tmp =
596 t_log_stretch_total(
i,
j) = t_log_u2_h1_tmp(
i,
j) / 2 + t_log_u(
i,
j);
599 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
608 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
610 t_h(
i,
j) = t_Omega(
i,
j) + t_u(
i,
j);
611 t_h_domega(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
616 t_levi_kirchoff_domega(
k,
l) = 0;
623 ++t_levi_kirchoff_dP;
624 ++t_levi_kirchoff_domega;
625 ++t_approx_P_adjont_dstretch;
626 ++t_approx_P_adjont_log_du;
627 ++t_approx_P_adjont_log_du_dP;
628 ++t_approx_P_adjont_log_du_domega;
639 ++t_log_stretch_total;
648 int nb_integration_pts = data.
getN().size1();
649 auto v = getVolume();
650 auto t_w = getFTensor0IntegrationWeight();
651 auto t_div_P = getFTensor1FromMat<3>(dataAtPts->divPAtPts);
652 auto t_s_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotAtPts);
653 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
654 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
655 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
656 dataAtPts->wL2DotDotAtPts.clear();
658 auto t_s_dot_dot_w = getFTensor1FromMat<3>(dataAtPts->wL2DotDotAtPts);
660 int nb_base_functions = data.
getN().size2();
664 auto get_ftensor1 = [](
auto &
v) {
669 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
671 auto t_nf = get_ftensor1(nF);
673 for (; bb != nb_dofs / 3; ++bb) {
674 t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
675 t_nf(
i) -=
a * t_row_base_fun * alphaW * t_s_dot_w(
i);
676 t_nf(
i) -=
a * t_row_base_fun * alphaRho * t_s_dot_dot_w(
i);
680 for (; bb != nb_base_functions; ++bb)
694 int nb_integration_pts = getGaussPts().size2();
695 auto v = getVolume();
696 auto t_w = getFTensor0IntegrationWeight();
697 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
698 int nb_base_functions = data.
getN().size2();
703 auto get_ftensor1 = [](
auto &
v) {
707 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
709 auto t_nf = get_ftensor1(nF);
711 for (; bb != nb_dofs / 3; ++bb) {
712 t_nf(
k) += (
a * t_row_base_fun) * t_levi_kirchoff(
k);
716 for (; bb != nb_base_functions; ++bb) {
728 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
729 CHKERR integratePiola(data);
732 CHKERR integratePolyconvexHencky(data);
734 CHKERR integrateHencky(data);
748 int nb_integration_pts = data.
getN().size1();
749 auto v = getVolume();
750 auto t_w = getFTensor0IntegrationWeight();
751 auto t_approx_P_adjont_log_du =
752 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
753 auto t_log_stretch_h1 =
754 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
756 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
758 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
764 auto get_ftensor2 = [](
auto &
v) {
766 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
769 int nb_base_functions = data.
getN().size2();
771 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
773 auto t_nf = get_ftensor2(nF);
777 t_D(
i,
j,
k,
l) * (t_log_stretch_h1(
k,
l) + alphaU * t_dot_log_u(
k,
l));
780 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
783 for (; bb != nb_dofs / 6; ++bb) {
784 t_nf(
L) += t_row_base_fun * t_residual(
L);
788 for (; bb != nb_base_functions; ++bb)
792 ++t_approx_P_adjont_log_du;
806 int nb_integration_pts = data.
getN().size1();
807 auto v = getVolume();
808 auto t_w = getFTensor0IntegrationWeight();
809 auto t_approx_P_adjont_log_du =
810 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
811 auto t_log_stretch_h1 =
812 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
814 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
816 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
822 auto get_ftensor2 = [](
auto &
v) {
824 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
827 constexpr
double nohat_k = 1. / 4;
828 constexpr
double hat_k = 1. / 8;
829 double mu = dataAtPts->mu;
830 double lambda = dataAtPts->lambda;
832 constexpr
double third = boost::math::constants::third<double>();
836 int nb_base_functions = data.
getN().size2();
838 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
840 auto t_nf = get_ftensor2(nF);
842 double log_det = t_log_stretch_h1(
i,
i);
843 double log_det2 = log_det * log_det;
846 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
849 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
850 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
853 A * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j))
861 alphaU * t_D(
i,
j,
k,
l) * t_log_stretch_h1(
k,
l);
865 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
868 for (; bb != nb_dofs / 6; ++bb) {
869 t_nf(
L) += t_row_base_fun * t_residual(
L);
873 for (; bb != nb_base_functions; ++bb)
877 ++t_approx_P_adjont_log_du;
891 int nb_integration_pts = data.
getN().size1();
892 auto v = getVolume();
893 auto t_w = getFTensor0IntegrationWeight();
894 auto t_approx_P_adjont_log_du =
895 getFTensor1FromMat<size_symm>(dataAtPts->adjointPdUAtPts);
896 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
898 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
900 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
906 auto get_ftensor2 = [](
auto &
v) {
908 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
911 int nb_base_functions = data.
getN().size2();
913 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
915 auto t_nf = get_ftensor2(nF);
918 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
928 a * (t_approx_P_adjont_log_du(
L) - t_Ldiff_u(
i,
j,
L) * t_P(
i,
j) -
929 t_L(
i,
j,
L) * t_viscous_P(
i,
j));
932 for (; bb != nb_dofs / 6; ++bb) {
933 t_nf(
L) += t_row_base_fun * t_residual(
L);
937 for (; bb != nb_base_functions; ++bb)
941 ++t_approx_P_adjont_log_du;
952 int nb_integration_pts = data.
getN().size1();
953 auto v = getVolume();
954 auto t_w = getFTensor0IntegrationWeight();
955 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
957 int nb_base_functions = data.
getN().size2() / 3;
964 auto get_ftensor1 = [](
auto &
v) {
969 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
971 auto t_nf = get_ftensor1(nF);
978 for (; bb != nb_dofs / 3; ++bb) {
979 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
984 for (; bb != nb_base_functions; ++bb)
997 int nb_integration_pts = data.
getN().size1();
998 auto v = getVolume();
999 auto t_w = getFTensor0IntegrationWeight();
1000 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
1002 int nb_base_functions = data.
getN().size2() / 9;
1009 auto get_ftensor0 = [](
auto &
v) {
1013 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1015 auto t_nf = get_ftensor0(nF);
1022 for (; bb != nb_dofs; ++bb) {
1023 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
1027 for (; bb != nb_base_functions; ++bb) {
1040 int nb_integration_pts = data.
getN().size1();
1041 auto v = getVolume();
1042 auto t_w = getFTensor0IntegrationWeight();
1043 auto t_w_l2 = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
1044 int nb_base_functions = data.
getN().size2() / 3;
1047 auto get_ftensor1 = [](
auto &
v) {
1052 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1054 auto t_nf = get_ftensor1(nF);
1056 for (; bb != nb_dofs / 3; ++bb) {
1057 double div_row_base = t_row_diff_base_fun(
i,
i);
1058 t_nf(
i) +=
a * div_row_base * t_w_l2(
i);
1060 ++t_row_diff_base_fun;
1062 for (; bb != nb_base_functions; ++bb) {
1063 ++t_row_diff_base_fun;
1077 for (
auto &bc : (*bcDispPtr)) {
1079 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1081 int nb_integration_pts = data.
getN().size1();
1082 auto t_normal = getFTensor1Normal();
1083 auto t_w = getFTensor0IntegrationWeight();
1084 int nb_base_functions = data.
getN().size2() / 3;
1089 auto get_ftensor1 = [](
auto &
v) {
1095 for (
auto &sm : scalingMethodsVec) {
1096 scale *= sm->getScale(getFEMethod()->ts_t);
1103 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1105 t_bc_res(
i) = t_bc_disp(
i);
1106 auto t_nf = get_ftensor1(nF);
1108 for (; bb != nb_dofs / 3; ++bb) {
1110 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_res(
i) * 0.5;
1114 for (; bb != nb_base_functions; ++bb)
1129 for (
auto &bc : (*bcRotPtr)) {
1131 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1133 int nb_integration_pts = data.
getN().size1();
1134 auto t_normal = getFTensor1Normal();
1135 auto t_w = getFTensor0IntegrationWeight();
1137 int nb_base_functions = data.
getN().size2() / 3;
1142 auto get_ftensor1 = [](
auto &
v) {
1150 double theta = bc.theta;
1151 theta *= getFEMethod()->ts_t;
1154 const double a = sqrt(t_normal(
i) * t_normal(
i));
1155 t_omega(
i) = t_normal(
i) * (theta /
a);
1157 auto t_coords = getFTensor1CoordsAtGaussPts();
1159 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1161 t_delta(
i) = t_center(
i) - t_coords(
i);
1163 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1165 auto t_nf = get_ftensor1(nF);
1167 for (; bb != nb_dofs / 3; ++bb) {
1168 t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1172 for (; bb != nb_base_functions; ++bb)
1187 int operator()(
int p_row,
int p_col,
int p_data)
const {
1188 return 2 * (p_data + 1);
1192 if (
ts_ctx == CTX_TSSETIFUNCTION) {
1202 CHKERR mField.loop_finite_elements(problemPtr->getName(),
"ESSENTIAL_BC",
1204 this->getCacheWeakPtr());
1216 auto t_normal = getFTensor1Normal();
1217 const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1219 t_unit_normal(
i) = t_normal(
i) / nrm2;
1221 int nb_integration_pts = data.
getN().size1();
1222 int nb_base_functions = data.
getN().size2() / 3;
1223 double ts_t = getFEMethod()->ts_t;
1225 auto integrate_matrix = [&]() {
1228 auto t_w = getFTensor0IntegrationWeight();
1229 matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
1233 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1236 for (; rr != nb_dofs / 3; ++rr) {
1237 const double a = t_row_base_fun(
i) * t_unit_normal(
i);
1239 for (
int cc = 0; cc != nb_dofs / 3; ++cc) {
1240 const double b = t_col_base_fun(
i) * t_unit_normal(
i);
1241 matPP(rr, cc) += t_w *
a * b;
1247 for (; rr != nb_base_functions; ++rr)
1256 auto integrate_rhs = [&](
auto &bc) {
1259 auto t_w = getFTensor0IntegrationWeight();
1260 vecPv.resize(3, nb_dofs / 3,
false);
1264 double ts_t = getFEMethod()->ts_t;
1266 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1268 for (; rr != nb_dofs / 3; ++rr) {
1269 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1270 for (
int dd = 0;
dd != 3; ++
dd)
1272 vecPv(
dd, rr) +=
t * bc.vals[
dd];
1275 for (; rr != nb_base_functions; ++rr)
1282 auto integrate_rhs_cook = [&](
auto &bc) {
1285 vecPv.resize(3, nb_dofs / 3,
false);
1288 auto t_w = getFTensor0IntegrationWeight();
1290 auto t_coords = getFTensor1CoordsAtGaussPts();
1292 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1294 auto calc_tau = [](
double y) {
1297 return -y * (y - 1) / 0.25;
1300 const double tau = calc_tau(t_coords(1));
1303 for (; rr != nb_dofs / 3; ++rr) {
1304 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1305 for (
int dd = 0;
dd != 3; ++
dd)
1307 vecPv(
dd, rr) +=
t * tau * bc.vals[
dd];
1311 for (; rr != nb_base_functions; ++rr)
1321 for (
auto &bc : *(bcFe.bcData)) {
1322 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1327 CHKERR integrate_matrix();
1328 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*")))
1329 CHKERR integrate_rhs_cook(bc);
1331 CHKERR integrate_rhs(bc);
1334 lapack_dposv(
'L', nb_dofs / 3, 3, &*matPP.data().begin(),
1335 nb_dofs / 3, &*vecPv.data().begin(), nb_dofs / 3);
1338 "The leading minor of order %i is "
1339 "not positive; definite;\nthe "
1340 "solution could not be computed",
1343 for (
int dd = 0;
dd != 3; ++
dd)
1345 for (
int rr = 0; rr != nb_dofs / 3; ++rr)
1357 int nb_integration_pts = row_data.
getN().size1();
1358 int row_nb_dofs = row_data.
getIndices().size();
1359 int col_nb_dofs = col_data.
getIndices().size();
1362 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2));
1365 auto v = getVolume();
1366 auto t_w = getFTensor0IntegrationWeight();
1367 int row_nb_base_functions = row_data.
getN().size2();
1369 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1372 for (; rr != row_nb_dofs / 3; ++rr) {
1374 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1375 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1376 double div_col_base = t_col_diff_base_fun(
i,
i);
1377 t_m(
i) +=
a * t_row_base_fun * div_col_base;
1379 ++t_col_diff_base_fun;
1383 for (; rr != row_nb_base_functions; ++rr)
1394 if (alphaW < std::numeric_limits<double>::epsilon() &&
1395 alphaRho < std::numeric_limits<double>::epsilon())
1398 const int nb_integration_pts = row_data.
getN().size1();
1399 const int row_nb_dofs = row_data.
getIndices().size();
1402 &
m(
r + 0,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 2,
c + 2)
1408 auto v = getVolume();
1409 auto t_w = getFTensor0IntegrationWeight();
1411 int row_nb_base_functions = row_data.
getN().size2();
1414 double ts_scale = alphaW * getTSa();
1415 if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon())
1416 ts_scale += alphaRho * getTSaa();
1418 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1419 double a =
v * t_w * ts_scale;
1422 for (; rr != row_nb_dofs / 3; ++rr) {
1425 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1426 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1427 const double b =
a * t_row_base_fun * t_col_base_fun;
1436 for (; rr != row_nb_base_functions; ++rr)
1451 int nb_integration_pts = row_data.
getN().size1();
1452 int row_nb_dofs = row_data.
getIndices().size();
1453 int col_nb_dofs = col_data.
getIndices().size();
1457 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
1459 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
1461 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
1463 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
1465 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
1467 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2));
1473 auto v = getVolume();
1474 auto t_w = getFTensor0IntegrationWeight();
1475 auto t_approx_P_adjont_log_du_dP =
1476 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1477 int row_nb_base_functions = row_data.
getN().size2();
1479 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1482 for (; rr != row_nb_dofs / 6; ++rr) {
1485 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1487 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1489 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
1497 for (; rr != row_nb_base_functions; ++rr)
1500 ++t_approx_P_adjont_log_du_dP;
1511 int nb_integration_pts = row_data.
getN().size1();
1512 int row_nb_dofs = row_data.
getIndices().size();
1513 int col_nb_dofs = col_data.
getIndices().size();
1516 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c), &
m(
r + 3,
c), &
m(
r + 4,
c),
1523 auto v = getVolume();
1524 auto t_w = getFTensor0IntegrationWeight();
1525 auto t_approx_P_adjont_log_du_dP =
1526 getFTensor3FromMat<3, 3, size_symm>(dataAtPts->adjointPdUdPAtPts);
1528 int row_nb_base_functions = row_data.
getN().size2();
1530 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1533 for (; rr != row_nb_dofs / 6; ++rr) {
1534 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1535 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1536 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1538 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
1545 for (; rr != row_nb_base_functions; ++rr)
1548 ++t_approx_P_adjont_log_du_dP;
1556 if (dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1557 CHKERR integratePiola(row_data, col_data);
1560 CHKERR integratePolyconvexHencky(row_data, col_data);
1562 CHKERR integrateHencky(row_data, col_data);
1577 int nb_integration_pts = row_data.
getN().size1();
1578 int row_nb_dofs = row_data.
getIndices().size();
1579 int col_nb_dofs = col_data.
getIndices().size();
1585 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
1586 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
1588 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
1589 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
1591 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
1592 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
1594 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
1595 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
1597 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
1598 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
1600 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
1601 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
1612 auto v = getVolume();
1613 auto t_w = getFTensor0IntegrationWeight();
1615 auto t_approx_P_adjont_dstretch =
1616 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1617 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1618 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1619 auto &nbUniq = dataAtPts->nbUniq;
1621 int row_nb_base_functions = row_data.
getN().size2();
1624 auto get_dP = [&]() {
1626 auto ts_a = getTSa();
1628 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1630 t_dP_tmp(
L,
J) = -(1 + alphaU * ts_a) *
1632 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
1634 if (EshelbianCore::stretchSelector >
LINEAR) {
1635 auto t_approx_P_adjont_dstretch =
1636 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1637 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1638 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1639 auto &nbUniq = dataAtPts->nbUniq;
1641 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1642 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1647 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
1648 t_approx_P_adjont_dstretch(
j,
i));
1653 t_dP(
L,
J) = t_L(
i,
j,
L) *
1654 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
1660 ++t_approx_P_adjont_dstretch;
1665 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1666 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1667 t_dP(
L,
J) = t_dP_tmp(
L,
J);
1672 return getFTensor2FromMat<size_symm, size_symm>(dP);
1675 auto t_dP = get_dP();
1676 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1680 for (; rr != row_nb_dofs / 6; ++rr) {
1682 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1683 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1684 const double b =
a * t_row_base_fun * t_col_base_fun;
1685 t_m(
L,
J) += b * t_dP(
L,
J);
1692 for (; rr != row_nb_base_functions; ++rr) {
1703 OpSpatialPhysical_du_du::integratePolyconvexHencky(
EntData &row_data,
1712 int nb_integration_pts = row_data.
getN().size1();
1713 int row_nb_dofs = row_data.
getIndices().size();
1714 int col_nb_dofs = col_data.
getIndices().size();
1720 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
1721 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
1723 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
1724 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
1726 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
1727 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
1729 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
1730 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
1732 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
1733 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
1735 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
1736 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
1747 auto v = getVolume();
1748 auto t_w = getFTensor0IntegrationWeight();
1750 int row_nb_base_functions = row_data.
getN().size2();
1753 auto get_dP = [&]() {
1755 auto ts_a = getTSa();
1757 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(dataAtPts->matD);
1759 constexpr
double nohat_k = 1. / 4;
1760 constexpr
double hat_k = 1. / 8;
1761 double mu = dataAtPts->mu;
1762 double lambda = dataAtPts->lambda;
1764 constexpr
double third = boost::math::constants::third<double>();
1768 auto t_approx_P_adjont_dstretch =
1769 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1770 auto t_log_stretch_h1 =
1771 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTotalTensorAtPts);
1772 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1773 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1774 auto &nbUniq = dataAtPts->nbUniq;
1776 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1777 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1779 double log_det = t_log_stretch_h1(
i,
i);
1780 double log_det2 = log_det * log_det;
1782 t_dev(
i,
j) = t_log_stretch_h1(
i,
j) -
t_kd(
i,
j) * (
third * log_det);
1783 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
1785 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
1786 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
1790 (
A * 2 * nohat_k) * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j));
1791 t_B_diff(
i,
j) = (B * 2 * hat_k) * log_det *
t_kd(
i,
j) +
1795 t_A_diff(
i,
j) * (t_dev(
m,
n) * t_diff_deviator(
m,
n,
k,
l))
1799 A * t_diff_deviator(
m,
n,
i,
j) * t_diff_deviator(
m,
n,
k,
l)
1805 t_dP(
L,
J) = -t_L(
i,
j,
L) *
1812 (alphaU * ts_a) * (t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)
1819 if (EshelbianCore::stretchSelector >
LINEAR) {
1821 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
1822 t_approx_P_adjont_dstretch(
j,
i));
1827 t_dP(
L,
J) += t_L(
i,
j,
L) *
1828 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
1834 ++t_approx_P_adjont_dstretch;
1840 return getFTensor2FromMat<size_symm, size_symm>(dP);
1843 auto t_dP = get_dP();
1844 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1848 for (; rr != row_nb_dofs / 6; ++rr) {
1850 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1851 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1852 const double b =
a * t_row_base_fun * t_col_base_fun;
1853 t_m(
L,
J) += b * t_dP(
L,
J);
1860 for (; rr != row_nb_base_functions; ++rr) {
1879 int nb_integration_pts = row_data.
getN().size1();
1880 int row_nb_dofs = row_data.
getIndices().size();
1881 int col_nb_dofs = col_data.
getIndices().size();
1887 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
1888 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
1890 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
1891 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
1893 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
1894 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
1896 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
1897 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
1899 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
1900 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
1902 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
1903 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
1914 auto v = getVolume();
1915 auto t_w = getFTensor0IntegrationWeight();
1917 auto get_dP = [&]() {
1919 auto ts_a = getTSa();
1921 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->PAtPts);
1922 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(dataAtPts->P_du);
1923 auto t_approx_P_adjont_dstretch =
1924 getFTensor2FromMat<3, 3>(dataAtPts->adjointPdstretchAtPts);
1926 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchDotTensorAtPts);
1927 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1929 getFTensor4DdgFromMat<3, 3, 1>(dataAtPts->diffStretchTensorAtPts);
1930 auto t_eigen_vals = getFTensor1FromMat<3>(dataAtPts->eigenVals);
1931 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(dataAtPts->eigenVecs);
1932 auto &nbUniq = dataAtPts->nbUniq;
1934 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(dP);
1935 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1938 t_deltaP(
i,
j) = t_approx_P_adjont_dstretch(
i,
j) - t_P(
i,
j);
1941 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
1943 -t_Ldiff_u(
i,
j,
L) * (r_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
1945 t_dP(
L,
J) -= (alphaU * ts_a) *
1946 (t_L(
i,
j,
L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
1948 if (EshelbianCore::stretchSelector >
LINEAR) {
1950 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
1951 t_deltaP_sym(
i,
j) /= 2.0;
1955 t_dP(
L,
J) += t_L(
i,
j,
L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
1961 ++t_approx_P_adjont_dstretch;
1969 return getFTensor2FromMat<size_symm, size_symm>(dP);
1972 int row_nb_base_functions = row_data.
getN().size2();
1975 auto t_dP = get_dP();
1976 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1980 for (; rr != row_nb_dofs / 6; ++rr) {
1982 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1983 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1984 const double b =
a * t_row_base_fun * t_col_base_fun;
1985 t_m(
L,
J) += b * t_dP(
L,
J);
1992 for (; rr != row_nb_base_functions; ++rr) {
2009 int row_nb_dofs = row_data.
getIndices().size();
2010 int col_nb_dofs = col_data.
getIndices().size();
2014 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2016 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2018 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2),
2020 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2),
2022 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2),
2024 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2)
2034 auto v = getVolume();
2035 auto t_w = getFTensor0IntegrationWeight();
2036 auto t_approx_P_adjont_log_du_domega =
2037 getFTensor2FromMat<3, size_symm>(dataAtPts->adjointPdUdOmegaAtPts);
2039 int row_nb_base_functions = row_data.
getN().size2();
2042 int nb_integration_pts = row_data.
getN().size1();
2044 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2048 for (; rr != row_nb_dofs / 6; ++rr) {
2050 auto t_m = get_ftensor3(
K, 6 * rr, 0);
2051 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2052 double v =
a * t_row_base_fun * t_col_base_fun;
2053 t_m(
L,
k) +=
v * t_approx_P_adjont_log_du_domega(
k,
L);
2060 for (; rr != row_nb_base_functions; ++rr)
2064 ++t_approx_P_adjont_log_du_domega;
2073 int nb_integration_pts = getGaussPts().size2();
2074 int row_nb_dofs = row_data.
getIndices().size();
2075 int col_nb_dofs = col_data.
getIndices().size();
2079 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2081 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2083 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
2093 auto v = getVolume();
2094 auto t_w = getFTensor0IntegrationWeight();
2095 auto t_levi_kirchoff_dP =
2096 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2098 int row_nb_base_functions = row_data.
getN().size2();
2101 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2104 for (; rr != row_nb_dofs / 3; ++rr) {
2105 double b =
a * t_row_base_fun;
2107 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2108 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2109 t_m(
k,
i) += b * (t_levi_kirchoff_dP(
i,
l,
k) * t_col_base_fun(
l));
2115 for (; rr != row_nb_base_functions; ++rr) {
2120 ++t_levi_kirchoff_dP;
2128 int nb_integration_pts = getGaussPts().size2();
2129 int row_nb_dofs = row_data.
getIndices().size();
2130 int col_nb_dofs = col_data.
getIndices().size();
2134 &
m(
r + 0,
c), &
m(
r + 1,
c), &
m(
r + 2,
c));
2143 auto v = getVolume();
2144 auto t_w = getFTensor0IntegrationWeight();
2145 auto t_levi_kirchoff_dP =
2146 getFTensor3FromMat<3, 3, 3>(dataAtPts->leviKirchhoffPAtPts);
2148 int row_nb_base_functions = row_data.
getN().size2();
2151 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2154 for (; rr != row_nb_dofs / 3; ++rr) {
2155 double b =
a * t_row_base_fun;
2156 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2157 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2158 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2159 t_m(
k) += b * (t_levi_kirchoff_dP(
i,
j,
k) * t_col_base_fun(
i,
j));
2166 for (; rr != row_nb_base_functions; ++rr) {
2170 ++t_levi_kirchoff_dP;
2178 int nb_integration_pts = getGaussPts().size2();
2179 int row_nb_dofs = row_data.
getIndices().size();
2180 int col_nb_dofs = col_data.
getIndices().size();
2184 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2186 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2188 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
2199 auto v = getVolume();
2200 auto t_w = getFTensor0IntegrationWeight();
2201 auto t_levi_kirchoff_domega =
2202 getFTensor2FromMat<3, 3>(dataAtPts->leviKirchhoffOmegaAtPts);
2203 int row_nb_base_functions = row_data.
getN().size2();
2205 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2208 for (; rr != row_nb_dofs / 3; ++rr) {
2209 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2210 const double b =
a * t_row_base_fun;
2212 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2213 t_m(
k,
l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(
k,
l);
2219 for (; rr != row_nb_base_functions; ++rr) {
2223 ++t_levi_kirchoff_domega;
2235 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2237 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2239 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
2244 int nb_integration_pts = getGaussPts().size2();
2245 int row_nb_dofs = row_data.
getIndices().size();
2246 int col_nb_dofs = col_data.
getIndices().size();
2248 auto v = getVolume();
2249 auto t_w = getFTensor0IntegrationWeight();
2250 int row_nb_base_functions = row_data.
getN().size2() / 3;
2260 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD);
2263 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2267 for (; rr != row_nb_dofs / 3; ++rr) {
2269 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2270 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2271 t_m(
i,
k) +=
a * t_row_base(
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
l));
2279 for (; rr != row_nb_base_functions; ++rr)
2289 OpSpatialConsistency_dBubble_dBubble::integrate(
EntData &row_data,
2293 int nb_integration_pts = getGaussPts().size2();
2294 int row_nb_dofs = row_data.
getIndices().size();
2295 int col_nb_dofs = col_data.
getIndices().size();
2297 auto v = getVolume();
2298 auto t_w = getFTensor0IntegrationWeight();
2299 int row_nb_base_functions = row_data.
getN().size2() / 9;
2309 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD);
2312 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2316 for (; rr != row_nb_dofs; ++rr) {
2318 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2320 a * (t_row_base(
i,
j) * (t_inv_D(
i,
j,
k,
l) * t_col_base(
k,
l)));
2327 for (; rr != row_nb_base_functions; ++rr)
2342 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2)
2347 int nb_integration_pts = getGaussPts().size2();
2348 int row_nb_dofs = row_data.
getIndices().size();
2349 int col_nb_dofs = col_data.
getIndices().size();
2351 auto v = getVolume();
2352 auto t_w = getFTensor0IntegrationWeight();
2353 int row_nb_base_functions = row_data.
getN().size2() / 9;
2363 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(dataAtPts->matInvD);
2366 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2369 auto t_m = get_ftensor1(
K, 0, 0);
2372 for (; rr != row_nb_dofs; ++rr) {
2374 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2375 t_m(
k) +=
a * (t_row_base(
i,
j) * t_inv_D(
i,
j,
k,
l)) * t_col_base(
l);
2383 for (; rr != row_nb_base_functions; ++rr)
2395 int nb_integration_pts = row_data.
getN().size1();
2396 int row_nb_dofs = row_data.
getIndices().size();
2397 int col_nb_dofs = col_data.
getIndices().size();
2402 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
2404 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
2406 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2)
2418 auto v = getVolume();
2419 auto t_w = getFTensor0IntegrationWeight();
2420 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2421 int row_nb_base_functions = row_data.
getN().size2() / 3;
2424 const double ts_a = getTSa();
2426 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2432 for (; rr != row_nb_dofs / 3; ++rr) {
2435 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2438 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2439 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2440 t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
2448 for (; rr != row_nb_base_functions; ++rr)
2457 OpSpatialConsistency_dBubble_domega::integrate(
EntData &row_data,
2461 int nb_integration_pts = row_data.
getN().size1();
2462 int row_nb_dofs = row_data.
getIndices().size();
2463 int col_nb_dofs = col_data.
getIndices().size();
2467 &
m(
r,
c + 0), &
m(
r,
c + 1), &
m(
r,
c + 2));
2477 auto v = getVolume();
2478 auto t_w = getFTensor0IntegrationWeight();
2479 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(dataAtPts->hdOmegaAtPts);
2480 int row_nb_base_functions = row_data.
getN().size2() / 9;
2483 const double ts_a = getTSa();
2485 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2489 for (; rr != row_nb_dofs; ++rr) {
2492 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
2495 auto t_m = get_ftensor2(
K, rr, 0);
2496 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2497 t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
2505 for (; rr != row_nb_base_functions; ++rr)
2518 auto create_tag = [
this](
const std::string tag_name,
const int size) {
2519 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2521 CHKERR postProcMesh.tag_get_handle(tag_name.c_str(), size, MB_TYPE_DOUBLE,
2522 th, MB_TAG_CREAT | MB_TAG_SPARSE,
2527 Tag th_w = create_tag(
"SpatialDisplacement", 3);
2528 Tag th_omega = create_tag(
"Omega", 3);
2529 Tag th_approxP = create_tag(
"Piola1Stress", 9);
2530 Tag th_sigma = create_tag(
"CauchyStress", 9);
2531 Tag th_log_u = create_tag(
"LogSpatialStretch", 9);
2532 Tag th_u = create_tag(
"SpatialStretch", 9);
2533 Tag th_h = create_tag(
"h", 9);
2534 Tag th_X = create_tag(
"X", 3);
2535 Tag th_detF = create_tag(
"detF", 1);
2536 Tag th_angular_momentum = create_tag(
"AngularMomentum", 3);
2538 Tag th_u_eig_vec = create_tag(
"SpatialStretchEigenVec", 9);
2539 Tag th_u_eig_vals = create_tag(
"SpatialStretchEigenVals", 3);
2540 Tag th_traction = create_tag(
"traction", 3);
2542 Tag th_disp = create_tag(
"H1Displacement", 3);
2543 Tag th_disp_error = create_tag(
"DisplacementError", 1);
2544 Tag th_lambda_disp = create_tag(
"ContactDisplacement", 3);
2546 auto t_w = getFTensor1FromMat<3>(dataAtPts->wL2AtPts);
2547 auto t_omega = getFTensor1FromMat<3>(dataAtPts->rotAxisAtPts);
2548 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2550 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
2551 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
2552 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
2553 auto t_approx_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2554 auto t_levi_kirchoff = getFTensor1FromMat<3>(dataAtPts->leviKirchhoffAtPts);
2555 auto t_coords = getFTensor1CoordsAtGaussPts();
2556 auto t_normal = getFTensor1NormalsAtGaussPts();
2557 auto t_disp = getFTensor1FromMat<3>(dataAtPts->wH1AtPts);
2558 auto t_lambda_disp = getFTensor1FromMat<3>(dataAtPts->contactL2AtPts);
2565 auto set_float_precision = [](
const double x) {
2566 if (std::abs(x) < std::numeric_limits<float>::epsilon())
2573 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
2575 v = set_float_precision(
v);
2576 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1, &
v);
2583 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2586 for (
auto &
a :
v.data())
2587 a = set_float_precision(
a);
2588 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
2589 &*
v.data().begin());
2597 &
m(0, 0), &
m(0, 1), &
m(0, 2),
2599 &
m(1, 0), &
m(1, 1), &
m(1, 2),
2601 &
m(2, 0), &
m(2, 1), &
m(2, 2));
2603 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2605 t_m(
i,
j) = t_d(
i,
j);
2606 for (
auto &
v :
m.data())
2607 v = set_float_precision(
v);
2608 CHKERR postProcMesh.tag_set_data(
th, &mapGaussPts[gg], 1,
2609 &*
m.data().begin());
2613 const auto nb_gauss_pts = getGaussPts().size2();
2614 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
2617 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
2620 CHKERR save_vec_tag(th_w, t_w, gg);
2621 CHKERR save_vec_tag(th_X, t_coords, gg);
2622 CHKERR save_vec_tag(th_omega, t_omega, gg);
2623 CHKERR save_vec_tag(th_traction, t_traction, gg);
2626 CHKERR save_mat_tag(th_h, t_h, gg);
2629 for (
int ii = 0; ii != 3; ++ii)
2630 for (
int jj = 0; jj != 3; ++jj)
2631 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2633 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2636 for (
int ii = 0; ii != 3; ++ii)
2637 for (
int jj = 0; jj != 3; ++jj)
2638 t_u_tmp(ii, jj) = t_u(ii, jj);
2640 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2641 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2642 CHKERR save_vec_tag(th_disp, t_disp, gg);
2643 CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
2645 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
2646 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2650 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
2651 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2652 CHKERR postProcMesh.tag_set_data(th_detF, &mapGaussPts[gg], 1, &jac);
2655 t_levi(
k) = t_levi_kirchoff(
k);
2656 CHKERR postProcMesh.tag_set_data(th_angular_momentum, &mapGaussPts[gg], 1,
2659 auto get_eiegn_vector_symmetric = [&](
auto &t_u) {
2662 for (
int ii = 0; ii != 3; ++ii)
2663 for (
int jj = 0; jj != 3; ++jj)
2664 t_m(ii, jj) = t_u(ii, jj);
2667 auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2670 CHKERR postProcMesh.tag_set_data(th_u_eig_vec, &mapGaussPts[gg], 1,
2671 &*
m.data().begin());
2672 CHKERR postProcMesh.tag_set_data(th_u_eig_vals, &mapGaussPts[gg], 1,
2673 &*eigen_values.data().begin());
2678 CHKERR get_eiegn_vector_symmetric(t_u);
2698 OpCalculateStrainEnergy::doWork(
int side, EntityType
type,
2701 if (
type == MBTET) {
2702 int nb_integration_pts = data.
getN().size1();
2703 auto v = getVolume();
2704 auto t_w = getFTensor0IntegrationWeight();
2705 auto t_P = getFTensor2FromMat<3, 3>(dataAtPts->approxPAtPts);
2706 auto t_h = getFTensor2FromMat<3, 3>(dataAtPts->hAtPts);
2711 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2712 const double a = t_w *
v;