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;
76 typedef typename std::remove_pointer<T>::type
Type;
79 typedef typename std::remove_pointer<T>::type
Type;
97 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
100 t_R(
i,
j) += t_Omega(
i,
j);
105 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
106 const auto a = sin(angle) / angle;
107 t_R(
i,
j) +=
a * t_Omega(
i,
j);
111 const auto ss_2 = sin(angle / 2.) / angle;
112 const auto b = 2. * ss_2 * ss_2;
113 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
132 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
134 t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
138 const auto ss = sin(angle);
139 const auto a = ss / angle;
140 t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
143 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
145 const auto angle2 = angle * angle;
146 const auto cc = cos(angle);
147 const auto diff_a = (angle * cc - ss) / angle2;
148 t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
152 const auto ss_2 = sin(angle / 2.);
153 const auto cc_2 = cos(angle / 2.);
154 const auto b = 2. * ss_2 * ss_2 / angle2;
156 b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
157 FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
159 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
161 diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
175 constexpr double eps = 1e-10;
176 for (
int l = 0;
l != 3; ++
l) {
178 t_omega_c(
i) = t_omega(
i);
179 t_omega_c(
l) += std::complex<double>(0,
eps);
181 for (
int i = 0;
i != 3; ++
i) {
182 for (
int j = 0;
j != 3; ++
j) {
183 for (
int k = 0;
k != 3; ++
k) {
184 t_diff2_R(
i,
j,
k,
l) = t_diff_R_c(
i,
j,
k).imag() /
eps;
194 static inline auto check(
const double &
a,
const double &b) {
195 constexpr double eps = std::numeric_limits<float>::epsilon();
203inline auto is_eq(
const double &
a,
const double &b) {
208 std::array<double, DIM> tmp;
209 std::copy(ptr, &ptr[DIM], tmp.begin());
210 std::sort(tmp.begin(), tmp.end());
211 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
212 constexpr double eps = std::numeric_limits<float>::epsilon();
214 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(),
is_eq));
222 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
224 int i = 0,
j = 1,
k = 2;
226 if (
is_eq(eig(0), eig(1))) {
230 }
else if (
is_eq(eig(0), eig(2))) {
234 }
else if (
is_eq(eig(1), eig(2))) {
241 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
243 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
245 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
253 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
265 int nb_integration_pts = getGaussPts().size2();
273 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
274 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
275 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
276 dataAtPts->leviKirchhoffAtPts.resize(3, nb_integration_pts,
false);
277 dataAtPts->leviKirchhoffPAtPts.resize(9 * 3, nb_integration_pts,
false);
278 dataAtPts->leviKirchhoffOmegaAtPts.resize(9, nb_integration_pts,
false);
280 dataAtPts->adjointPdstretchAtPts.resize(9, nb_integration_pts,
false);
285 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
286 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
287 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
288 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
289 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
290 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
291 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
293 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
295 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
296 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
298 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->hdLogStretchAtPts);
299 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchhoffAtPts);
300 auto t_levi_kirchoff_dP =
301 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchhoffPAtPts);
302 auto t_levi_kirchoff_domega =
303 getFTensor2FromMat<3, 3>(
dataAtPts->leviKirchhoffOmegaAtPts);
304 auto t_approx_P_adjont_dstretch =
305 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
306 auto t_approx_P_adjont_log_du =
307 getFTensor1FromMat<size_symm>(
dataAtPts->adjointPdUAtPts);
308 auto t_approx_P_adjont_log_du_dP =
309 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjointPdUdPAtPts);
310 auto t_approx_P_adjont_log_du_domega =
311 getFTensor2FromMat<3, size_symm>(
dataAtPts->adjointPdUdOmegaAtPts);
312 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
313 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
314 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(
dataAtPts->diffRotMatAtPts);
316 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
317 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
318 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
320 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
321 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
322 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
325 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
326 auto t_log_stretch_total =
327 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
330 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
341 auto calulate_rotation = [&]() {
345 t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
346 t_R(
i,
j) = t0(
i,
j);
349 auto calulate_streach = [&]() {
350 eigen_vec(
i,
j) = t_log_u(
i,
j);
353 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
354 if (nbUniq[gg] < 3) {
355 sort_eigen_vals<3>(eig, eigen_vec);
357 t_eigen_vals(
i) = eig(
i);
358 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
361 t_u(
i,
j) = t_u_tmp(
i,
j);
365 t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
366 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
376 t_Ru(
i,
m) = t_R(
i,
l) * t_u(
l,
m);
377 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
378 t_h_domega(
i,
j,
k) = (t_diff_R(
i,
l,
k) * t_u(
l,
m)) * t_h1(
m,
j);
379 t_h_dlog_u(
i,
j,
L) = (t_R(
i,
l) * t_Ldiff_u(
l,
m,
L)) * t_h1(
m,
j);
382 t_approx_P_intermidiate(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
383 t_approx_P_adjont_dstretch(
l,
m) =
384 t_approx_P_intermidiate(
i,
m) * t_R(
i,
l);
387 levi_civita(
l,
m,
k) * (t_approx_P_adjont_dstretch(
l,
m));
388 t_levi_kirchoff_dP(
i,
j,
k) =
389 (levi_civita(
l,
m,
k) * t_R(
i,
l)) * t_h1(
m,
j);
390 t_levi_kirchoff_domega(
k,
n) =
391 levi_civita(
l,
m,
k) *
392 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
394 t_approx_P_adjont_log_du(
L) =
395 t_Ldiff_u(
l,
m,
L) * t_approx_P_adjont_dstretch(
l,
m);
396 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
397 t_approx_P_adjont_log_du_domega(
n,
L) =
399 (t_approx_P_intermidiate(
i,
m) * t_diff_R(
i,
l,
n));
408 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
409 t_Ru(
i,
m) = t_Omega(
i,
m) + t_u(
i,
m);
410 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
411 t_h_domega(
i,
j,
k) = 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);
419 levi_civita(
i,
m,
k) * (t_approx_P_adjont_dstretch(
i,
m));
420 t_levi_kirchoff_dP(
i,
j,
k) = levi_civita(
i,
m,
k) * t_h1(
m,
j);
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);
444 t_levi_kirchoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
445 t_levi_kirchoff_dP(
i,
j,
k) = levi_civita(
i,
j,
k);
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);
467 for (
int ii = 0; ii != 3; ++ii) {
468 eig(ii) = std::max(eig(ii),
470 std::numeric_limits<double>::min_exponent)));
475 auto t_log_u2_h1_tmp =
481 t_log_stretch_total(
i,
j) = t_log_u2_h1_tmp(
i,
j) / 2 + t_log_u(
i,
j);
484 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
492 ++t_levi_kirchoff_dP;
493 ++t_levi_kirchoff_domega;
494 ++t_approx_P_adjont_dstretch;
495 ++t_approx_P_adjont_log_du;
496 ++t_approx_P_adjont_log_du_dP;
497 ++t_approx_P_adjont_log_du_domega;
508 ++t_log_stretch_total;
517 int nb_integration_pts = data.
getN().size1();
518 auto v = getVolume();
519 auto t_w = getFTensor0IntegrationWeight();
520 auto t_div_P = getFTensor1FromMat<3>(
dataAtPts->divPAtPts);
521 auto t_s_dot_w = getFTensor1FromMat<3>(
dataAtPts->wL2DotAtPts);
522 if (
dataAtPts->wL2DotDotAtPts.size1() == 0 &&
523 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
524 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
527 auto t_s_dot_dot_w = getFTensor1FromMat<3>(
dataAtPts->wL2DotDotAtPts);
529 int nb_base_functions = data.
getN().size2();
533 auto get_ftensor1 = [](
auto &
v) {
538 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
540 auto t_nf = get_ftensor1(
nF);
542 for (; bb != nb_dofs / 3; ++bb) {
543 t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
544 t_nf(
i) -=
a * t_row_base_fun *
alphaW * t_s_dot_w(
i);
545 t_nf(
i) -=
a * t_row_base_fun *
alphaRho * t_s_dot_dot_w(
i);
549 for (; bb != nb_base_functions; ++bb)
563 int nb_integration_pts = getGaussPts().size2();
564 auto v = getVolume();
565 auto t_w = getFTensor0IntegrationWeight();
566 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchhoffAtPts);
567 int nb_base_functions = data.
getN().size2();
572 auto get_ftensor1 = [](
auto &
v) {
576 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
578 auto t_nf = get_ftensor1(
nF);
580 for (; bb != nb_dofs / 3; ++bb) {
581 t_nf(
k) += (
a * t_row_base_fun) * t_levi_kirchoff(
k);
585 for (; bb != nb_base_functions; ++bb) {
597 if (
dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
617 int nb_integration_pts = data.
getN().size1();
618 auto v = getVolume();
619 auto t_w = getFTensor0IntegrationWeight();
620 auto t_approx_P_adjont_log_du =
621 getFTensor1FromMat<size_symm>(
dataAtPts->adjointPdUAtPts);
622 auto t_log_streach_h1 =
623 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
625 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
627 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
633 auto get_ftensor2 = [](
auto &
v) {
635 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
638 int nb_base_functions = data.
getN().size2();
640 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
642 auto t_nf = get_ftensor2(
nF);
646 t_D(
i,
j,
k,
l) * (t_log_streach_h1(
k,
l) +
alphaU * t_dot_log_u(
k,
l));
649 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
652 for (; bb != nb_dofs / 6; ++bb) {
653 t_nf(
L) += t_row_base_fun * t_residual(
L);
657 for (; bb != nb_base_functions; ++bb)
661 ++t_approx_P_adjont_log_du;
675 int nb_integration_pts = data.
getN().size1();
676 auto v = getVolume();
677 auto t_w = getFTensor0IntegrationWeight();
678 auto t_approx_P_adjont_log_du =
679 getFTensor1FromMat<size_symm>(
dataAtPts->adjointPdUAtPts);
680 auto t_log_streach_h1 =
681 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
683 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
685 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
691 auto get_ftensor2 = [](
auto &
v) {
693 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
696 constexpr double nohat_k = 1. / 4;
697 constexpr double hat_k = 1. / 8;
701 constexpr double third = boost::math::constants::third<double>();
705 int nb_base_functions = data.
getN().size2();
707 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
709 auto t_nf = get_ftensor2(
nF);
711 double log_det = t_log_streach_h1(
i,
i);
712 double log_det2 = log_det * log_det;
715 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
718 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
719 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
722 A * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j))
734 a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
737 for (; bb != nb_dofs / 6; ++bb) {
738 t_nf(
L) += t_row_base_fun * t_residual(
L);
742 for (; bb != nb_base_functions; ++bb)
746 ++t_approx_P_adjont_log_du;
760 int nb_integration_pts = data.
getN().size1();
761 auto v = getVolume();
762 auto t_w = getFTensor0IntegrationWeight();
763 auto t_approx_P_adjont_log_du =
764 getFTensor1FromMat<size_symm>(
dataAtPts->adjointPdUAtPts);
765 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
767 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
769 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
775 auto get_ftensor2 = [](
auto &
v) {
777 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
780 int nb_base_functions = data.
getN().size2();
782 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
784 auto t_nf = get_ftensor2(
nF);
787 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
797 a * (t_approx_P_adjont_log_du(
L) - t_Ldiff_u(
i,
j,
L) * t_P(
i,
j) -
798 t_L(
i,
j,
L) * t_viscous_P(
i,
j));
801 for (; bb != nb_dofs / 6; ++bb) {
802 t_nf(
L) += t_row_base_fun * t_residual(
L);
806 for (; bb != nb_base_functions; ++bb)
810 ++t_approx_P_adjont_log_du;
821 int nb_integration_pts = data.
getN().size1();
822 auto v = getVolume();
823 auto t_w = getFTensor0IntegrationWeight();
824 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
825 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
827 int nb_base_functions = data.
getN().size2() / 3;
834 auto get_ftensor1 = [](
auto &
v) {
839 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
841 auto t_nf = get_ftensor1(
nF);
848 for (; bb != nb_dofs / 3; ++bb) {
849 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
854 for (; bb != nb_base_functions; ++bb)
868 int nb_integration_pts = data.
getN().size1();
869 auto v = getVolume();
870 auto t_w = getFTensor0IntegrationWeight();
871 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
872 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
874 int nb_base_functions = data.
getN().size2() / 9;
881 auto get_ftensor0 = [](
auto &
v) {
885 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
887 auto t_nf = get_ftensor0(
nF);
894 for (; bb != nb_dofs; ++bb) {
895 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
899 for (; bb != nb_base_functions; ++bb) {
913 int nb_integration_pts = data.
getN().size1();
914 auto v = getVolume();
915 auto t_w = getFTensor0IntegrationWeight();
916 auto t_w_l2 = getFTensor1FromMat<3>(
dataAtPts->wL2AtPts);
917 int nb_base_functions = data.
getN().size2() / 3;
920 auto get_ftensor1 = [](
auto &
v) {
925 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
927 auto t_nf = get_ftensor1(
nF);
929 for (; bb != nb_dofs / 3; ++bb) {
930 double div_row_base = t_row_diff_base_fun(
i,
i);
931 t_nf(
i) +=
a * div_row_base * t_w_l2(
i);
933 ++t_row_diff_base_fun;
935 for (; bb != nb_base_functions; ++bb) {
936 ++t_row_diff_base_fun;
950 for (
auto &bc : (*bcDispPtr)) {
952 if (bc.faces.find(fe_ent) != bc.faces.end()) {
954 int nb_integration_pts = data.
getN().size1();
955 auto t_normal = getFTensor1Normal();
956 auto t_w = getFTensor0IntegrationWeight();
957 int nb_base_functions = data.
getN().size2() / 3;
962 auto get_ftensor1 = [](
auto &
v) {
969 scale *= sm->getScale(getFEMethod()->ts_t);
976 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
978 t_bc_res(
i) = t_bc_disp(
i);
979 auto t_nf = get_ftensor1(
nF);
981 for (; bb != nb_dofs / 3; ++bb) {
983 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_res(
i) * 0.5;
987 for (; bb != nb_base_functions; ++bb)
1002 for (
auto &bc : (*bcRotPtr)) {
1004 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1006 int nb_integration_pts = data.
getN().size1();
1007 auto t_normal = getFTensor1Normal();
1008 auto t_w = getFTensor0IntegrationWeight();
1010 int nb_base_functions = data.
getN().size2() / 3;
1015 auto get_ftensor1 = [](
auto &
v) {
1023 double theta = bc.theta;
1024 theta *= getFEMethod()->ts_t;
1027 const double a = sqrt(t_normal(
i) * t_normal(
i));
1028 t_omega(
i) = t_normal(
i) * (theta /
a);
1030 auto t_coords = getFTensor1CoordsAtGaussPts();
1032 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1034 t_delta(
i) = t_center(
i) - t_coords(
i);
1036 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
1038 auto t_nf = get_ftensor1(
nF);
1040 for (; bb != nb_dofs / 3; ++bb) {
1041 t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
1045 for (; bb != nb_base_functions; ++bb)
1060 int operator()(
int p_row,
int p_col,
int p_data)
const {
1061 return 2 * (p_data + 1);
1065 if (
ts_ctx == CTX_TSSETIFUNCTION) {
1089 auto t_normal = getFTensor1Normal();
1090 const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
1092 t_unit_normal(
i) = t_normal(
i) / nrm2;
1094 int nb_integration_pts = data.
getN().size1();
1095 int nb_base_functions = data.
getN().size2() / 3;
1096 double ts_t = getFEMethod()->ts_t;
1098 auto integrate_matrix = [&]() {
1101 auto t_w = getFTensor0IntegrationWeight();
1102 matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
1106 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1109 for (; rr != nb_dofs / 3; ++rr) {
1110 const double a = t_row_base_fun(
i) * t_unit_normal(
i);
1112 for (
int cc = 0; cc != nb_dofs / 3; ++cc) {
1113 const double b = t_col_base_fun(
i) * t_unit_normal(
i);
1114 matPP(rr, cc) += t_w *
a * b;
1120 for (; rr != nb_base_functions; ++rr)
1129 auto integrate_rhs = [&](
auto &bc) {
1132 auto t_w = getFTensor0IntegrationWeight();
1133 vecPv.resize(3, nb_dofs / 3,
false);
1137 double ts_t = getFEMethod()->ts_t;
1139 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1141 for (; rr != nb_dofs / 3; ++rr) {
1142 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1143 for (
int dd = 0; dd != 3; ++dd)
1145 vecPv(dd, rr) +=
t * bc.vals[dd];
1148 for (; rr != nb_base_functions; ++rr)
1155 auto integrate_rhs_cook = [&](
auto &bc) {
1158 vecPv.resize(3, nb_dofs / 3,
false);
1161 auto t_w = getFTensor0IntegrationWeight();
1163 auto t_coords = getFTensor1CoordsAtGaussPts();
1165 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1167 auto calc_tau = [](
double y) {
1170 return -y * (y - 1) / 0.25;
1173 const double tau = calc_tau(t_coords(1));
1176 for (; rr != nb_dofs / 3; ++rr) {
1177 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1178 for (
int dd = 0; dd != 3; ++dd)
1180 vecPv(dd, rr) +=
t * tau * bc.vals[dd];
1184 for (; rr != nb_base_functions; ++rr)
1195 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1200 CHKERR integrate_matrix();
1201 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*")))
1202 CHKERR integrate_rhs_cook(bc);
1204 CHKERR integrate_rhs(bc);
1208 nb_dofs / 3, &*
vecPv.data().begin(), nb_dofs / 3);
1211 "The leading minor of order %i is "
1212 "not positive; definite;\nthe "
1213 "solution could not be computed",
1216 for (
int dd = 0; dd != 3; ++dd)
1218 for (
int rr = 0; rr != nb_dofs / 3; ++rr)
1230 int nb_integration_pts = row_data.
getN().size1();
1231 int row_nb_dofs = row_data.
getIndices().size();
1232 int col_nb_dofs = col_data.
getIndices().size();
1233 auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
1235 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
1238 auto v = getVolume();
1239 auto t_w = getFTensor0IntegrationWeight();
1240 int row_nb_base_functions = row_data.
getN().size2();
1242 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1245 for (; rr != row_nb_dofs / 3; ++rr) {
1247 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1248 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1249 double div_col_base = t_col_diff_base_fun(
i,
i);
1250 t_m(
i) +=
a * t_row_base_fun * div_col_base;
1252 ++t_col_diff_base_fun;
1256 for (; rr != row_nb_base_functions; ++rr)
1267 if (
alphaW < std::numeric_limits<double>::epsilon() &&
1268 alphaRho < std::numeric_limits<double>::epsilon())
1271 const int nb_integration_pts = row_data.
getN().size1();
1272 const int row_nb_dofs = row_data.
getIndices().size();
1273 auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
1275 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
1281 auto v = getVolume();
1282 auto t_w = getFTensor0IntegrationWeight();
1284 int row_nb_base_functions = row_data.
getN().size2();
1287 double ts_scale =
alphaW * getTSa();
1288 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon())
1291 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1292 double a =
v * t_w * ts_scale;
1295 for (; rr != row_nb_dofs / 3; ++rr) {
1298 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1299 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1300 const double b =
a * t_row_base_fun * t_col_base_fun;
1309 for (; rr != row_nb_base_functions; ++rr)
1324 int nb_integration_pts = row_data.
getN().size1();
1325 int row_nb_dofs = row_data.
getIndices().size();
1326 int col_nb_dofs = col_data.
getIndices().size();
1327 auto get_ftensor3 = [](MatrixDouble &
m,
const int r,
const int c) {
1330 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1332 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1334 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
1336 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
1338 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
1340 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
1346 auto v = getVolume();
1347 auto t_w = getFTensor0IntegrationWeight();
1348 auto t_approx_P_adjont_log_du_dP =
1349 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjointPdUdPAtPts);
1350 int row_nb_base_functions = row_data.
getN().size2();
1352 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1355 for (; rr != row_nb_dofs / 6; ++rr) {
1358 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1360 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1362 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
1370 for (; rr != row_nb_base_functions; ++rr)
1373 ++t_approx_P_adjont_log_du_dP;
1384 int nb_integration_pts = row_data.
getN().size1();
1385 int row_nb_dofs = row_data.
getIndices().size();
1386 int col_nb_dofs = col_data.
getIndices().size();
1387 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
1389 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
1396 auto v = getVolume();
1397 auto t_w = getFTensor0IntegrationWeight();
1398 auto t_approx_P_adjont_log_du_dP =
1399 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjointPdUdPAtPts);
1401 int row_nb_base_functions = row_data.
getN().size2();
1403 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1406 for (; rr != row_nb_dofs / 6; ++rr) {
1407 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1408 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1409 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1411 a * (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
1418 for (; rr != row_nb_base_functions; ++rr)
1421 ++t_approx_P_adjont_log_du_dP;
1429 if (
dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1450 int nb_integration_pts = row_data.
getN().size1();
1451 int row_nb_dofs = row_data.
getIndices().size();
1452 int col_nb_dofs = col_data.
getIndices().size();
1454 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
1458 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1459 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1461 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1462 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1464 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1465 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1467 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1468 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1470 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1471 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1473 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1474 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1485 auto v = getVolume();
1486 auto t_w = getFTensor0IntegrationWeight();
1488 auto t_approx_P_adjont_dstretch =
1489 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
1490 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1491 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1494 int row_nb_base_functions = row_data.
getN().size2();
1497 auto get_dP = [&]() {
1499 auto ts_a = getTSa();
1501 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
1503 t_dP_tmp(
L,
J) = -(1 +
alphaU * ts_a) *
1505 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
1509 auto t_approx_P_adjont_dstretch =
1510 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
1511 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1512 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1515 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(
dP);
1516 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1521 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
1522 t_approx_P_adjont_dstretch(
j,
i));
1527 t_dP(
L,
J) = t_L(
i,
j,
L) *
1528 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
1534 ++t_approx_P_adjont_dstretch;
1539 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(
dP);
1540 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1541 t_dP(
L,
J) = t_dP_tmp(
L,
J);
1546 return getFTensor2FromMat<size_symm, size_symm>(
dP);
1549 auto t_dP = get_dP();
1550 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1554 for (; rr != row_nb_dofs / 6; ++rr) {
1556 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1557 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1558 const double b =
a * t_row_base_fun * t_col_base_fun;
1559 t_m(
L,
J) += b * t_dP(
L,
J);
1566 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();
1590 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
1594 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1595 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1597 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1598 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1600 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1601 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1603 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1604 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1606 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1607 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1609 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1610 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1621 auto v = getVolume();
1622 auto t_w = getFTensor0IntegrationWeight();
1624 int row_nb_base_functions = row_data.
getN().size2();
1627 auto get_dP = [&]() {
1629 auto ts_a = getTSa();
1631 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
1633 constexpr double nohat_k = 1. / 4;
1634 constexpr double hat_k = 1. / 8;
1638 constexpr double third = boost::math::constants::third<double>();
1642 auto t_approx_P_adjont_dstretch =
1643 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
1644 auto t_log_streach_h1 =
1645 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
1646 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1647 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1650 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(
dP);
1651 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1653 double log_det = t_log_streach_h1(
i,
i);
1654 double log_det2 = log_det * log_det;
1656 t_dev(
i,
j) = t_log_streach_h1(
i,
j) -
t_kd(
i,
j) * (
third * log_det);
1657 double dev_norm2 = t_dev(
i,
j) * t_dev(
i,
j);
1659 auto A = 2 *
mu * std::exp(nohat_k * dev_norm2);
1660 auto B =
lambda * std::exp(hat_k * log_det2) * log_det;
1664 (A * 2 * nohat_k) * (t_dev(
k,
l) * t_diff_deviator(
k,
l,
i,
j));
1665 t_B_diff(
i,
j) = (
B * 2 * hat_k) * log_det *
t_kd(
i,
j) +
1669 t_A_diff(
i,
j) * (t_dev(
m,
n) * t_diff_deviator(
m,
n,
k,
l))
1673 A * t_diff_deviator(
m,
n,
i,
j) * t_diff_deviator(
m,
n,
k,
l)
1679 t_dP(
L,
J) = -t_L(
i,
j,
L) *
1695 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
1696 t_approx_P_adjont_dstretch(
j,
i));
1701 t_dP(
L,
J) += t_L(
i,
j,
L) *
1702 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) *
1708 ++t_approx_P_adjont_dstretch;
1714 return getFTensor2FromMat<size_symm, size_symm>(
dP);
1717 auto t_dP = get_dP();
1718 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1722 for (; rr != row_nb_dofs / 6; ++rr) {
1724 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1725 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1726 const double b =
a * t_row_base_fun * t_col_base_fun;
1727 t_m(
L,
J) += b * t_dP(
L,
J);
1734 for (; rr != row_nb_base_functions; ++rr) {
1753 int nb_integration_pts = row_data.
getN().size1();
1754 int row_nb_dofs = row_data.
getIndices().size();
1755 int col_nb_dofs = col_data.
getIndices().size();
1757 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
1761 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1762 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1764 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1765 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1767 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1768 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1770 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1771 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1773 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1774 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1776 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1777 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1788 auto v = getVolume();
1789 auto t_w = getFTensor0IntegrationWeight();
1791 auto get_dP = [&]() {
1793 auto ts_a = getTSa();
1795 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
1796 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(
dataAtPts->P_du);
1797 auto t_approx_P_adjont_dstretch =
1798 getFTensor2FromMat<3, 3>(
dataAtPts->adjointPdstretchAtPts);
1800 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
1801 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
1803 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
1804 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1805 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1808 auto t_dP = getFTensor2FromMat<size_symm, size_symm>(
dP);
1809 for (
auto gg = 0; gg != nb_integration_pts; ++gg) {
1812 t_deltaP(
i,
j) = t_approx_P_adjont_dstretch(
i,
j) - t_P(
i,
j);
1815 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
1817 -t_Ldiff_u(
i,
j,
L) * (r_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
1820 (t_L(
i,
j,
L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
1824 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
1825 t_deltaP_sym(
i,
j) /= 2.0;
1829 t_dP(
L,
J) += t_L(
i,
j,
L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
1835 ++t_approx_P_adjont_dstretch;
1843 return getFTensor2FromMat<size_symm, size_symm>(
dP);
1846 int row_nb_base_functions = row_data.
getN().size2();
1849 auto t_dP = get_dP();
1850 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1854 for (; rr != row_nb_dofs / 6; ++rr) {
1856 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1857 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1858 const double b =
a * t_row_base_fun * t_col_base_fun;
1859 t_m(
L,
J) += b * t_dP(
L,
J);
1866 for (; rr != row_nb_base_functions; ++rr) {
1883 int row_nb_dofs = row_data.
getIndices().size();
1884 int col_nb_dofs = col_data.
getIndices().size();
1885 auto get_ftensor3 = [](MatrixDouble &
m,
const int r,
const int c) {
1888 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1890 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1892 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
1894 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
1896 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
1898 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
1908 auto v = getVolume();
1909 auto t_w = getFTensor0IntegrationWeight();
1910 auto t_approx_P_adjont_log_du_domega =
1911 getFTensor2FromMat<3, size_symm>(
dataAtPts->adjointPdUdOmegaAtPts);
1913 int row_nb_base_functions = row_data.
getN().size2();
1916 int nb_integration_pts = row_data.
getN().size1();
1918 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1922 for (; rr != row_nb_dofs / 6; ++rr) {
1924 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1925 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1926 double v =
a * t_row_base_fun * t_col_base_fun;
1927 t_m(
L,
k) +=
v * t_approx_P_adjont_log_du_domega(
k,
L);
1934 for (; rr != row_nb_base_functions; ++rr)
1938 ++t_approx_P_adjont_log_du_domega;
1947 int nb_integration_pts = getGaussPts().size2();
1948 int row_nb_dofs = row_data.
getIndices().size();
1949 int col_nb_dofs = col_data.
getIndices().size();
1950 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
1953 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1955 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1957 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1967 auto v = getVolume();
1968 auto t_w = getFTensor0IntegrationWeight();
1969 auto t_levi_kirchoff_dP =
1970 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchhoffPAtPts);
1972 int row_nb_base_functions = row_data.
getN().size2();
1975 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1978 for (; rr != row_nb_dofs / 3; ++rr) {
1979 double b =
a * t_row_base_fun;
1981 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1982 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1983 t_m(
k,
i) += b * (t_levi_kirchoff_dP(
i,
l,
k) * t_col_base_fun(
l));
1989 for (; rr != row_nb_base_functions; ++rr) {
1994 ++t_levi_kirchoff_dP;
2002 int nb_integration_pts = getGaussPts().size2();
2003 int row_nb_dofs = row_data.
getIndices().size();
2004 int col_nb_dofs = col_data.
getIndices().size();
2006 auto get_ftensor1 = [](MatrixDouble &
m,
const int r,
const int c) {
2008 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
2017 auto v = getVolume();
2018 auto t_w = getFTensor0IntegrationWeight();
2019 auto t_levi_kirchoff_dP =
2020 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchhoffPAtPts);
2022 int row_nb_base_functions = row_data.
getN().size2();
2025 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2028 for (; rr != row_nb_dofs / 3; ++rr) {
2029 double b =
a * t_row_base_fun;
2030 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
2031 auto t_m = get_ftensor1(
K, 3 * rr, 0);
2032 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
2033 t_m(
k) += b * (t_levi_kirchoff_dP(
i,
j,
k) * t_col_base_fun(
i,
j));
2040 for (; rr != row_nb_base_functions; ++rr) {
2044 ++t_levi_kirchoff_dP;
2052 int nb_integration_pts = getGaussPts().size2();
2053 int row_nb_dofs = row_data.
getIndices().size();
2054 int col_nb_dofs = col_data.
getIndices().size();
2055 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
2058 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2060 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2062 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2073 auto v = getVolume();
2074 auto t_w = getFTensor0IntegrationWeight();
2075 auto t_levi_kirchoff_domega =
2076 getFTensor2FromMat<3, 3>(
dataAtPts->leviKirchhoffOmegaAtPts);
2077 int row_nb_base_functions = row_data.
getN().size2();
2079 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2082 for (; rr != row_nb_dofs / 3; ++rr) {
2083 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2084 const double b =
a * t_row_base_fun;
2086 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2087 t_m(
k,
l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(
k,
l);
2093 for (; rr != row_nb_base_functions; ++rr) {
2097 ++t_levi_kirchoff_domega;
2106 int nb_integration_pts = row_data.
getN().size1();
2107 int row_nb_dofs = row_data.
getIndices().size();
2108 int col_nb_dofs = col_data.
getIndices().size();
2110 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
2113 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
2115 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
2117 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
2129 auto v = getVolume();
2130 auto t_w = getFTensor0IntegrationWeight();
2131 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
2132 int row_nb_base_functions = row_data.
getN().size2() / 3;
2135 const double ts_a = getTSa();
2137 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2143 for (; rr != row_nb_dofs / 3; ++rr) {
2146 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
2149 auto t_m = get_ftensor2(
K, 3 * rr, 0);
2150 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2151 t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
2159 for (; rr != row_nb_base_functions; ++rr)
2172 int nb_integration_pts = row_data.
getN().size1();
2173 int row_nb_dofs = row_data.
getIndices().size();
2174 int col_nb_dofs = col_data.
getIndices().size();
2176 auto get_ftensor2 = [](MatrixDouble &
m,
const int r,
const int c) {
2178 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
2188 auto v = getVolume();
2189 auto t_w = getFTensor0IntegrationWeight();
2190 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
2191 int row_nb_base_functions = row_data.
getN().size2() / 9;
2194 const double ts_a = getTSa();
2196 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2200 for (; rr != row_nb_dofs; ++rr) {
2203 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
2206 auto t_m = get_ftensor2(
K, rr, 0);
2207 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2208 t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
2216 for (; rr != row_nb_base_functions; ++rr)
2229 auto create_tag = [
this](
const std::string tag_name,
const int size) {
2230 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
2233 th, MB_TAG_CREAT | MB_TAG_SPARSE,
2238 Tag th_w = create_tag(
"SpatialDisplacement", 3);
2239 Tag th_omega = create_tag(
"Omega", 3);
2240 Tag th_approxP = create_tag(
"Piola1Stress", 9);
2241 Tag th_sigma = create_tag(
"CauchyStress", 9);
2242 Tag th_log_u = create_tag(
"LogSpatialStretch", 9);
2243 Tag th_u = create_tag(
"SpatialStretch", 9);
2244 Tag th_h = create_tag(
"h", 9);
2245 Tag th_X = create_tag(
"X", 3);
2246 Tag th_detF = create_tag(
"detF", 1);
2247 Tag th_angular_momentum = create_tag(
"AngularMomentum", 3);
2249 Tag th_u_eig_vec = create_tag(
"SpatialStretchEigenVec", 9);
2250 Tag th_u_eig_vals = create_tag(
"SpatialStretchEigenVals", 3);
2251 Tag th_traction = create_tag(
"traction", 3);
2253 Tag th_disp = create_tag(
"H1Displacement", 3);
2254 Tag th_disp_error = create_tag(
"DisplacementError", 1);
2255 Tag th_lambda_disp = create_tag(
"ContactDisplacement", 3);
2257 auto t_w = getFTensor1FromMat<3>(
dataAtPts->wL2AtPts);
2258 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
2259 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
2261 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
2262 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
2263 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
2264 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
2265 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchhoffAtPts);
2266 auto t_coords = getFTensor1CoordsAtGaussPts();
2267 auto t_normal = getFTensor1NormalsAtGaussPts();
2268 auto t_disp = getFTensor1FromMat<3>(
dataAtPts->wH1AtPts);
2269 auto t_lambda_disp = getFTensor1FromMat<3>(
dataAtPts->contactL2AtPts);
2276 auto set_float_precision = [](
const double x) {
2277 if (std::abs(x) < std::numeric_limits<float>::epsilon())
2284 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
2286 v = set_float_precision(
v);
2294 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2297 for (
auto &
a :
v.data())
2298 a = set_float_precision(
a);
2300 &*
v.data().begin());
2306 MatrixDouble3by3
m(3, 3);
2308 &
m(0, 0), &
m(0, 1), &
m(0, 2),
2310 &
m(1, 0), &
m(1, 1), &
m(1, 2),
2312 &
m(2, 0), &
m(2, 1), &
m(2, 2));
2314 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
2316 t_m(
i,
j) = t_d(
i,
j);
2317 for (
auto &
v :
m.data())
2318 v = set_float_precision(
v);
2320 &*
m.data().begin());
2324 const auto nb_gauss_pts = getGaussPts().size2();
2325 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
2328 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
2331 CHKERR save_vec_tag(th_w, t_w, gg);
2332 CHKERR save_vec_tag(th_X, t_coords, gg);
2333 CHKERR save_vec_tag(th_omega, t_omega, gg);
2334 CHKERR save_vec_tag(th_traction, t_traction, gg);
2337 CHKERR save_mat_tag(th_h, t_h, gg);
2340 for (
int ii = 0; ii != 3; ++ii)
2341 for (
int jj = 0; jj != 3; ++jj)
2342 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2344 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2347 for (
int ii = 0; ii != 3; ++ii)
2348 for (
int jj = 0; jj != 3; ++jj)
2349 t_u_tmp(ii, jj) = t_u(ii, jj);
2351 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2352 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2353 CHKERR save_vec_tag(th_disp, t_disp, gg);
2354 CHKERR save_vec_tag(th_lambda_disp, t_lambda_disp, gg);
2356 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
2357 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2361 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
2362 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2366 t_levi(
k) = t_levi_kirchoff(
k);
2370 auto get_eiegn_vector_symmetric = [&](
auto &t_u) {
2373 for (
int ii = 0; ii != 3; ++ii)
2374 for (
int jj = 0; jj != 3; ++jj)
2375 t_m(ii, jj) = t_u(ii, jj);
2377 VectorDouble3 eigen_values(3);
2378 auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2382 &*
m.data().begin());
2384 &*eigen_values.data().begin());
2389 CHKERR get_eiegn_vector_symmetric(t_u);
2412 if (type == MBTET) {
2413 int nb_integration_pts = data.
getN().size1();
2414 auto v = getVolume();
2415 auto t_w = getFTensor0IntegrationWeight();
2416 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
2417 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
2422 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2423 const double a = t_w *
v;
2425 (*energy) +=
a * t_P(
i,
J) * t_h(
i,
J);
Eshelbian plasticity interface.
Kronecker Delta class symmetric.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
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)
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
auto get_uniq_nb(double *ptr)
auto diff_deviator(FTensor::Ddg< double, 3, 3 > &&t_diff_stress)
auto get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
static constexpr auto size_symm
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
auto get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
auto is_eq(const double &a, const double &b)
auto get_diff2_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega, RotSelector rotSelector=LARGE_ROT)
Tensors class implemented by Walter Landry.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
static boost::function< double(const double)> dd_f
static enum RotSelector gradApperoximator
static boost::function< double(const double)> inv_f
static enum StretchSelector stretchSelector
static boost::function< double(const double)> d_f
static enum RotSelector rotSelector
static boost::function< double(const double)> f
MoFEMErrorCode preProcess()
boost::shared_ptr< TractionBcVec > bcData
MoFEM::Interface & mField
MatrixDouble K
local tangent matrix
VectorDouble nF
local right hand side vector
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 doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &data)
std::vector< boost::shared_ptr< ScalingMethod > > scalingMethodsVec
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePolyconvexHencky(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(EntData &data)
MoFEMErrorCode integratePolyconvexHencky(EntData &data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode integrateHencky(EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
static auto check(const double &a, const double &b)
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
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::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.