12#include <boost/math/constants/constants.hpp>
49 typedef typename std::remove_pointer<T>::type
Type;
52 typedef typename std::remove_pointer<T>::type
Type;
70 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
73 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
75 t_R(
i,
j) += t_Omega(
i,
j);
79 const auto a = sin(angle) / angle;
80 t_R(
i,
j) +=
a * t_Omega(
i,
j);
84 const auto ss_2 = sin(angle / 2.);
85 const auto b = 2. * ss_2 * ss_2 / (angle * angle);
86 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
105 sqrt(t_omega(
i) * t_omega(
i)) + std::numeric_limits<double>::epsilon();
107 t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
j,
k);
111 const auto ss = sin(angle);
112 const auto a = ss / angle;
113 t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
116 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
118 const auto angle2 = angle * angle;
119 const auto cc = cos(angle);
120 const auto diff_a = (angle * cc - ss) / angle2;
121 t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
125 const auto ss_2 = sin(angle / 2.);
126 const auto cc_2 = cos(angle / 2.);
127 const auto b = 2. * ss_2 * ss_2 / angle2;
129 b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
130 FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
132 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
134 diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
148 constexpr double eps = 1e-10;
149 for (
int l = 0;
l != 3; ++
l) {
151 t_omega_c(
i) = t_omega(
i);
152 t_omega_c(
l) += std::complex<double>(0,
eps);
154 for (
int i = 0;
i != 3; ++
i) {
155 for (
int j = 0;
j != 3; ++
j) {
156 for (
int k = 0;
k != 3; ++
k) {
157 t_diff2_R(
i,
j,
k,
l) = t_diff_R_c(
i,
j,
k).imag() /
eps;
167 static inline auto check(
const double &a,
const double &b) {
168 constexpr double eps = std::numeric_limits<float>::epsilon();
176inline auto is_eq(
const double &a,
const double &b) {
181 std::array<double, DIM> tmp;
182 std::copy(ptr, &ptr[DIM], tmp.begin());
183 std::sort(tmp.begin(), tmp.end());
184 isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
185 constexpr double eps = std::numeric_limits<float>::epsilon();
187 return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(),
is_eq));
195 std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
197 int i = 0,
j = 1,
k = 2;
199 if (
is_eq(eig(0), eig(1))) {
203 }
else if (
is_eq(eig(0), eig(2))) {
207 }
else if (
is_eq(eig(1), eig(2))) {
214 eigen_vec(
i, 0), eigen_vec(
i, 1), eigen_vec(
i, 2),
216 eigen_vec(
j, 0), eigen_vec(
j, 1), eigen_vec(
j, 2),
218 eigen_vec(
k, 0), eigen_vec(
k, 1), eigen_vec(
k, 2)};
226 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
238 int nb_integration_pts = getGaussPts().size2();
246 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
247 dataAtPts->hdOmegaAtPts.resize(9 * 3, nb_integration_pts,
false);
248 dataAtPts->hdLogStretchAtPts.resize(9 * 6, nb_integration_pts,
false);
249 dataAtPts->leviKirchoffAtPts.resize(3, nb_integration_pts,
false);
250 dataAtPts->leviKirchoffdPAtPts.resize(9 * 3, nb_integration_pts,
false);
251 dataAtPts->leviKirchoffdOmegaAtPts.resize(9, nb_integration_pts,
false);
253 dataAtPts->adjontPdstretchAtPts.resize(9, nb_integration_pts,
false);
258 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
259 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
260 dataAtPts->stretchTensorAtPts.resize(6, nb_integration_pts,
false);
261 dataAtPts->diffStretchTensorAtPts.resize(36, nb_integration_pts,
false);
262 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
263 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
264 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
266 dataAtPts->logStretchTotalTensorAtPts.resize(6, nb_integration_pts,
false);
268 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
269 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
270 auto t_h_dlog_u = getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->hdLogStretchAtPts);
271 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchoffAtPts);
272 auto t_levi_kirchoff_dP =
273 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchoffdPAtPts);
274 auto t_levi_kirchoff_domega =
275 getFTensor2FromMat<3, 3>(
dataAtPts->leviKirchoffdOmegaAtPts);
276 auto t_approx_P_adjont_dstretch =
277 getFTensor2FromMat<3, 3>(
dataAtPts->adjontPdstretchAtPts);
278 auto t_approx_P_adjont_log_du =
279 getFTensor1FromMat<size_symm>(
dataAtPts->adjontPdUAtPts);
280 auto t_approx_P_adjont_log_du_dP =
281 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjontPdUdPAtPts);
282 auto t_approx_P_adjont_log_du_domega =
283 getFTensor2FromMat<3, size_symm>(
dataAtPts->adjontPdUdOmegaAtPts);
284 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
285 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
286 auto t_diff_R = getFTensor3FromMat<3, 3, 3>(
dataAtPts->diffRotMatAtPts);
288 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
289 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
290 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
292 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
293 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
294 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
297 auto t_grad_h1 = getFTensor2FromMat<3, 3>(
dataAtPts->wGradH1AtPts);
298 auto t_log_stretch_total =
299 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
302 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
313 auto calulate_rotation = [&]() {
317 t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
318 t_R(
i,
j) = t0(
i,
j);
321 auto calulate_streach = [&]() {
322 eigen_vec(
i,
j) = t_log_u(
i,
j);
325 nbUniq[gg] = get_uniq_nb<3>(&eig(0));
326 if (nbUniq[gg] < 3) {
327 sort_eigen_vals<3>(eig, eigen_vec);
329 t_eigen_vals(
i) = eig(
i);
330 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
333 t_u(
i,
j) = t_u_tmp(
i,
j);
337 t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
338 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
347 t_Ru(
i,
m) = t_R(
i,
l) * t_u(
l,
m);
348 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
349 t_h_domega(
i,
j,
k) = (t_diff_R(
i,
l,
k) * t_u(
l,
m)) * t_h1(
m,
j);
350 t_h_dlog_u(
i,
j,
L) = (t_R(
i,
l) * t_Ldiff_u(
l,
m,
L)) * t_h1(
m,
j);
352 t_approx_P_intermidiata(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
353 t_approx_P_adjont_dstretch(
l,
m) =
354 t_approx_P_intermidiata(
i,
m) * t_R(
i,
l);
357 levi_civita(
l,
m,
k) * (t_approx_P_adjont_dstretch(
l,
m));
358 t_levi_kirchoff_dP(
i,
j,
k) =
359 (levi_civita(
l,
m,
k) * t_R(
i,
l)) * t_h1(
m,
j);
360 t_levi_kirchoff_domega(
k,
n) =
361 levi_civita(
l,
m,
k) *
362 (t_approx_P_intermidiata(
i,
m) * t_diff_R(
i,
l,
n));
364 t_approx_P_adjont_log_du(
L) =
365 t_Ldiff_u(
l,
m,
L) * t_approx_P_adjont_dstretch(
l,
m);
366 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
367 t_approx_P_adjont_log_du_domega(
n,
L) =
369 (t_approx_P_intermidiata(
i,
m) * t_diff_R(
i,
l,
n));
376 t_h(
i,
j) = t_Ru(
i,
m) * t_h1(
m,
j);
377 t_h_domega(
i,
j,
k) = t_diff_R(
i,
m,
k) * t_h1(
m,
j);
378 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
m,
L) * t_h1(
m,
j);
380 t_approx_P_intermidiata(
i,
m) = t_approx_P(
i,
j) * t_h1(
m,
j);
381 t_approx_P_adjont_dstretch(
i,
m) = t_approx_P_intermidiata(
i,
m);
384 levi_civita(
i,
m,
k) * (t_approx_P_adjont_dstretch(
i,
m));
385 t_levi_kirchoff_dP(
i,
j,
k) = levi_civita(
i,
m,
k) * t_h1(
m,
j);
386 t_levi_kirchoff_domega(
k,
n) = 0;
388 t_approx_P_adjont_log_du(
L) =
389 t_Ldiff_u(
i,
m,
L) * t_approx_P_adjont_dstretch(
i,
m);
390 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
391 t_approx_P_adjont_log_du_domega(
n,
L) = 0;
399 t_h_domega(
i,
j,
k) = t_diff_R(
i,
j,
k);
400 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
402 t_levi_kirchoff(
k) = levi_civita(
i,
j,
k) * t_approx_P(
i,
j);
403 t_levi_kirchoff_dP(
i,
j,
k) = levi_civita(
i,
j,
k);
404 t_levi_kirchoff_domega(
k,
l) = 0;
406 t_approx_P_adjont_dstretch(
i,
j) = t_approx_P(
i,
j);
407 t_approx_P_adjont_log_du(
L) =
408 t_Ldiff_u(
i,
j,
L) * t_approx_P_adjont_dstretch(
i,
j);
409 t_approx_P_adjont_log_du_dP(
i,
j,
L) = t_h_dlog_u(
i,
j,
L);
410 t_approx_P_adjont_log_du_domega(
m,
L) = 0;
416 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
417 eigen_vec(
i,
j) = t_C_h1(
i,
j);
423 for (
int ii = 0; ii != 3; ++ii) {
424 eig(ii) = std::max(eig(ii),
426 std::numeric_limits<double>::min_exponent)));
431 auto t_log_u2_h1_tmp =
437 t_log_stretch_total(
i,
j) = t_log_u2_h1_tmp(
i,
j) / 2 + t_log_u(
i,
j);
440 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
448 ++t_levi_kirchoff_dP;
449 ++t_levi_kirchoff_domega;
450 ++t_approx_P_adjont_dstretch;
451 ++t_approx_P_adjont_log_du;
452 ++t_approx_P_adjont_log_du_dP;
453 ++t_approx_P_adjont_log_du_domega;
464 ++t_log_stretch_total;
473 int nb_integration_pts = data.
getN().size1();
474 auto v = getVolume();
475 auto t_w = getFTensor0IntegrationWeight();
476 auto t_div_P = getFTensor1FromMat<3>(
dataAtPts->divPAtPts);
477 auto t_s_dot_w = getFTensor1FromMat<3>(
dataAtPts->wL2DotAtPts);
478 if (
dataAtPts->wL2DotDotAtPts.size1() == 0 &&
479 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
480 dataAtPts->wL2DotDotAtPts.resize(3, nb_integration_pts);
483 auto t_s_dot_dot_w = getFTensor1FromMat<3>(
dataAtPts->wL2DotDotAtPts);
485 int nb_base_functions = data.
getN().size2();
489 auto get_ftensor1 = [](
auto &
v) {
494 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
496 auto t_nf = get_ftensor1(
nF);
498 for (; bb != nb_dofs / 3; ++bb) {
499 t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
500 t_nf(
i) -=
a * t_row_base_fun *
alphaW * t_s_dot_w(
i);
501 t_nf(
i) -=
a * t_row_base_fun *
alphaRho * t_s_dot_dot_w(
i);
505 for (; bb != nb_base_functions; ++bb)
519 int nb_integration_pts = getGaussPts().size2();
520 auto v = getVolume();
521 auto t_w = getFTensor0IntegrationWeight();
522 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchoffAtPts);
523 int nb_base_functions = data.
getN().size2();
528 auto get_ftensor1 = [](
auto &
v) {
532 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
534 auto t_nf = get_ftensor1(
nF);
536 for (; bb != nb_dofs / 3; ++bb) {
537 t_nf(
k) += (
a * t_row_base_fun) * t_levi_kirchoff(
k);
541 for (; bb != nb_base_functions; ++bb) {
553 if (
dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
569 int nb_integration_pts = data.
getN().size1();
570 auto v = getVolume();
571 auto t_w = getFTensor0IntegrationWeight();
572 auto t_approx_P_adjont_log_du =
573 getFTensor1FromMat<size_symm>(
dataAtPts->adjontPdUAtPts);
574 auto t_log_streach_h1 =
575 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTotalTensorAtPts);
577 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
579 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
585 auto get_ftensor2 = [](
auto &
v) {
587 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
590 int nb_base_functions = data.
getN().size2();
592 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
594 auto t_nf = get_ftensor2(
nF);
598 t_D(
i,
j,
k,
l) * (t_log_streach_h1(
k,
l) +
alphaU * t_dot_log_u(
k,
l));
600 t_residual(
L) =
a * (t_approx_P_adjont_log_du(
L) - t_L(
i,
j,
L) * t_T(
i,
j));
603 for (; bb != nb_dofs / 6; ++bb) {
604 t_nf(
L) += t_row_base_fun * t_residual(
L);
608 for (; bb != nb_base_functions; ++bb)
612 ++t_approx_P_adjont_log_du;
626 int nb_integration_pts = data.
getN().size1();
627 auto v = getVolume();
628 auto t_w = getFTensor0IntegrationWeight();
629 auto t_approx_P_adjont_log_du =
630 getFTensor1FromMat<size_symm>(
dataAtPts->adjontPdUAtPts);
631 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
633 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
635 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
641 auto get_ftensor2 = [](
auto &
v) {
643 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
646 int nb_base_functions = data.
getN().size2();
648 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
650 auto t_nf = get_ftensor2(
nF);
653 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
k,
l) * t_L(
k,
l,
L);
663 a * (t_approx_P_adjont_log_du(
L) - t_Ldiff_u(
i,
j,
L) * t_P(
i,
j) -
664 t_L(
i,
j,
L) * t_viscous_P(
i,
j));
667 for (; bb != nb_dofs / 6; ++bb) {
668 t_nf(
L) += t_row_base_fun * t_residual(
L);
672 for (; bb != nb_base_functions; ++bb)
676 ++t_approx_P_adjont_log_du;
687 int nb_integration_pts = data.
getN().size1();
688 auto v = getVolume();
689 auto t_w = getFTensor0IntegrationWeight();
690 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
691 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
693 int nb_base_functions = data.
getN().size2() / 3;
700 auto get_ftensor1 = [](
auto &
v) {
705 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
707 auto t_nf = get_ftensor1(
nF);
714 for (; bb != nb_dofs / 3; ++bb) {
715 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
720 for (; bb != nb_base_functions; ++bb)
734 int nb_integration_pts = data.
getN().size1();
735 auto v = getVolume();
736 auto t_w = getFTensor0IntegrationWeight();
737 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
738 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
740 int nb_base_functions = data.
getN().size2() / 9;
747 auto get_ftensor0 = [](
auto &
v) {
751 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
753 auto t_nf = get_ftensor0(
nF);
760 for (; bb != nb_dofs; ++bb) {
761 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
765 for (; bb != nb_base_functions; ++bb) {
779 int nb_integration_pts = data.
getN().size1();
780 auto v = getVolume();
781 auto t_w = getFTensor0IntegrationWeight();
782 auto t_w_l2 = getFTensor1FromMat<3>(
dataAtPts->wL2AtPts);
783 int nb_base_functions = data.
getN().size2() / 3;
786 auto get_ftensor1 = [](
auto &
v) {
791 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
793 auto t_nf = get_ftensor1(
nF);
795 for (; bb != nb_dofs / 3; ++bb) {
796 double div_row_base = t_row_diff_base_fun(
i,
i);
797 t_nf(
i) +=
a * div_row_base * t_w_l2(
i);
799 ++t_row_diff_base_fun;
801 for (; bb != nb_base_functions; ++bb) {
802 ++t_row_diff_base_fun;
818 if (bc.faces.find(fe_ent) != bc.faces.end()) {
820 int nb_integration_pts = data.
getN().size1();
821 auto t_normal = getFTensor1Normal();
822 auto t_w = getFTensor0IntegrationWeight();
823 int nb_base_functions = data.
getN().size2() / 3;
828 auto get_ftensor1 = [](
auto &
v) {
835 t_bc_disp(
i) *= getFEMethod()->ts_t;
837 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
839 t_bc_res(
i) = t_bc_disp(
i);
840 auto t_nf = get_ftensor1(
nF);
842 for (; bb != nb_dofs / 3; ++bb) {
844 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_res(
i) * 0.5;
848 for (; bb != nb_base_functions; ++bb)
865 if (bc.faces.find(fe_ent) != bc.faces.end()) {
867 int nb_integration_pts = data.
getN().size1();
868 auto t_normal = getFTensor1Normal();
869 auto t_w = getFTensor0IntegrationWeight();
871 int nb_base_functions = data.
getN().size2() / 3;
876 auto get_ftensor1 = [](
auto &
v) {
884 double theta = bc.theta;
885 theta *= getFEMethod()->ts_t;
888 const double a = sqrt(t_normal(
i) * t_normal(
i));
889 t_omega(
i) = t_normal(
i) * (theta /
a);
891 auto t_coords = getFTensor1CoordsAtGaussPts();
893 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
895 t_delta(
i) = t_center(
i) - t_coords(
i);
897 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
899 auto t_nf = get_ftensor1(
nF);
901 for (; bb != nb_dofs / 3; ++bb) {
902 t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
906 for (; bb != nb_base_functions; ++bb)
921 int operator()(
int p_row,
int p_col,
int p_data)
const {
922 return 2 * (p_data + 1);
926 if (
ts_ctx == CTX_TSSETIFUNCTION) {
938 this->getCacheWeakPtr());
950 auto t_normal = getFTensor1Normal();
951 const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
953 t_unit_normal(
i) = t_normal(
i) / nrm2;
955 int nb_integration_pts = data.
getN().size1();
956 int nb_base_functions = data.
getN().size2() / 3;
957 double ts_t = getFEMethod()->ts_t;
959 auto integrate_matrix = [&]() {
962 auto t_w = getFTensor0IntegrationWeight();
963 matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
967 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
970 for (; rr != nb_dofs / 3; ++rr) {
971 const double a = t_row_base_fun(
i) * t_unit_normal(
i);
973 for (
int cc = 0; cc != nb_dofs / 3; ++cc) {
974 const double b = t_col_base_fun(
i) * t_unit_normal(
i);
975 matPP(rr, cc) += t_w *
a * b;
981 for (; rr != nb_base_functions; ++rr)
990 auto integrate_rhs = [&](
auto &bc) {
993 auto t_w = getFTensor0IntegrationWeight();
994 vecPv.resize(3, nb_dofs / 3,
false);
998 double ts_t = getFEMethod()->ts_t;
1000 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1002 for (; rr != nb_dofs / 3; ++rr) {
1003 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1004 for (
int dd = 0; dd != 3; ++dd)
1006 vecPv(dd, rr) +=
t * bc.vals[dd];
1009 for (; rr != nb_base_functions; ++rr)
1016 auto integrate_rhs_cook = [&](
auto &bc) {
1019 vecPv.resize(3, nb_dofs / 3,
false);
1022 auto t_w = getFTensor0IntegrationWeight();
1024 auto t_coords = getFTensor1CoordsAtGaussPts();
1026 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1028 auto calc_tau = [](
double y) {
1031 return -y * (y - 1) / 0.25;
1034 const double tau = calc_tau(t_coords(1));
1037 for (; rr != nb_dofs / 3; ++rr) {
1038 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1039 for (
int dd = 0; dd != 3; ++dd)
1041 vecPv(dd, rr) +=
t * tau * bc.vals[dd];
1045 for (; rr != nb_base_functions; ++rr)
1056 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1061 CHKERR integrate_matrix();
1062 if (std::regex_match(bc.blockName, std::regex(
".*COOK.*")))
1063 CHKERR integrate_rhs_cook(bc);
1065 CHKERR integrate_rhs(bc);
1069 nb_dofs / 3, &*
vecPv.data().begin(), nb_dofs / 3);
1072 "The leading minor of order %i is "
1073 "not positive; definite;\nthe "
1074 "solution could not be computed",
1077 for (
int dd = 0; dd != 3; ++dd)
1079 for (
int rr = 0; rr != nb_dofs / 3; ++rr)
1091 int nb_integration_pts = row_data.
getN().size1();
1092 int row_nb_dofs = row_data.
getIndices().size();
1093 int col_nb_dofs = col_data.
getIndices().size();
1094 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
1096 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
1099 auto v = getVolume();
1100 auto t_w = getFTensor0IntegrationWeight();
1101 int row_nb_base_functions = row_data.
getN().size2();
1103 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1106 for (; rr != row_nb_dofs / 3; ++rr) {
1108 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1109 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1110 double div_col_base = t_col_diff_base_fun(
i,
i);
1111 t_m(
i) +=
a * t_row_base_fun * div_col_base;
1113 ++t_col_diff_base_fun;
1117 for (; rr != row_nb_base_functions; ++rr)
1128 if (
alphaW < std::numeric_limits<double>::epsilon() &&
1129 alphaRho < std::numeric_limits<double>::epsilon())
1132 const int nb_integration_pts = row_data.
getN().size1();
1133 const int row_nb_dofs = row_data.
getIndices().size();
1134 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
1136 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
1142 auto v = getVolume();
1143 auto t_w = getFTensor0IntegrationWeight();
1145 int row_nb_base_functions = row_data.
getN().size2();
1148 double ts_scale =
alphaW * getTSa();
1149 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon())
1152 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1153 double a =
v * t_w * ts_scale;
1156 for (; rr != row_nb_dofs / 3; ++rr) {
1159 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1160 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
1161 const double b =
a * t_row_base_fun * t_col_base_fun;
1170 for (; rr != row_nb_base_functions; ++rr)
1185 int nb_integration_pts = row_data.
getN().size1();
1186 int row_nb_dofs = row_data.
getIndices().size();
1187 int col_nb_dofs = col_data.
getIndices().size();
1188 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
1191 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1193 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1195 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
1197 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
1199 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
1201 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
1207 auto v = getVolume();
1208 auto t_w = getFTensor0IntegrationWeight();
1209 auto t_approx_P_adjont_log_du_dP =
1210 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjontPdUdPAtPts);
1211 int row_nb_base_functions = row_data.
getN().size2();
1213 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1216 for (; rr != row_nb_dofs / 6; ++rr) {
1219 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1221 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1223 (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
j)) *
1231 for (; rr != row_nb_base_functions; ++rr)
1234 ++t_approx_P_adjont_log_du_dP;
1245 int nb_integration_pts = row_data.
getN().size1();
1246 int row_nb_dofs = row_data.
getIndices().size();
1247 int col_nb_dofs = col_data.
getIndices().size();
1248 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1250 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
1257 auto v = getVolume();
1258 auto t_w = getFTensor0IntegrationWeight();
1259 auto t_approx_P_adjont_log_du_dP =
1260 getFTensor3FromMat<3, 3, size_symm>(
dataAtPts->adjontPdUdPAtPts);
1262 int row_nb_base_functions = row_data.
getN().size2();
1264 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1267 for (; rr != row_nb_dofs / 6; ++rr) {
1268 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1269 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1270 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1272 (t_approx_P_adjont_log_du_dP(
i,
j,
L) * t_col_base_fun(
i,
j)) *
1279 for (; rr != row_nb_base_functions; ++rr)
1282 ++t_approx_P_adjont_log_du_dP;
1290 if (
dataAtPts->physicsPtr->dependentVariablesPiola.size()) {
1307 int nb_integration_pts = row_data.
getN().size1();
1308 int row_nb_dofs = row_data.
getIndices().size();
1309 int col_nb_dofs = col_data.
getIndices().size();
1311 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1315 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1316 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1318 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1319 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1321 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1322 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1324 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1325 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1327 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1328 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1330 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1331 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1346 auto v = getVolume();
1347 auto t_w = getFTensor0IntegrationWeight();
1349 auto t_approx_P_adjont_dstretch =
1350 getFTensor2FromMat<3, 3>(
dataAtPts->adjontPdstretchAtPts);
1352 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
1353 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
1355 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
1356 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1357 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1360 int row_nb_base_functions = row_data.
getN().size2();
1363 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(
dataAtPts->matD);
1365 const double ts_a = getTSa();
1366 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1375 t_sym(
i,
j) = (t_approx_P_adjont_dstretch(
i,
j) ||
1376 t_approx_P_adjont_dstretch(
j,
i));
1383 ((t_diff2_uP2(
i,
j,
k,
l) + t_diff2_uP2(
k,
l,
i,
j)) * t_L(
k,
l,
J)) /
1389 t_dP(
L,
J) -= (1 +
alphaU * ts_a) *
1391 ((t_D(
i,
j,
m,
n) * t_diff(
m,
n,
k,
l)) * t_L(
k,
l,
J)));
1394 for (; rr != row_nb_dofs / 6; ++rr) {
1396 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1397 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1398 const double b =
a * t_row_base_fun * t_col_base_fun;
1399 t_m(
L,
J) += b * t_dP(
L,
J);
1406 for (; rr != row_nb_base_functions; ++rr) {
1411 ++t_approx_P_adjont_dstretch;
1430 int nb_integration_pts = row_data.
getN().size1();
1431 int row_nb_dofs = row_data.
getIndices().size();
1432 int col_nb_dofs = col_data.
getIndices().size();
1434 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1438 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1439 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1441 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1442 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1444 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1445 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1447 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1448 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1450 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1451 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1453 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1454 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1465 auto v = getVolume();
1466 auto t_w = getFTensor0IntegrationWeight();
1467 auto r_P_du = getFTensor4FromMat<3, 3, 3, 3>(
dataAtPts->P_du);
1469 auto t_approx_P_adjont_dstretch =
1470 getFTensor2FromMat<3, 3>(
dataAtPts->adjontPdstretchAtPts);
1471 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
1473 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchDotTensorAtPts);
1474 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
1476 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStretchTensorAtPts);
1477 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1478 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1481 int row_nb_base_functions = row_data.
getN().size2();
1484 const double ts_a = getTSa();
1485 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1489 t_deltaP(
i,
j) = t_approx_P_adjont_dstretch(
i,
j) - t_P(
i,
j);
1497 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
1498 t_deltaP_sym(
i,
j) /= 2.0;
1502 t_dP(
L,
J) = t_L(
i,
j,
L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
1508 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
1510 t_Ldiff_u(
i,
j,
L) * (r_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
1516 for (; rr != row_nb_dofs / 6; ++rr) {
1518 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1519 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1520 const double b =
a * t_row_base_fun * t_col_base_fun;
1521 t_m(
L,
J) += b * t_dP(
L,
J);
1528 for (; rr != row_nb_base_functions; ++rr) {
1535 ++t_approx_P_adjont_dstretch;
1552 int row_nb_dofs = row_data.
getIndices().size();
1553 int col_nb_dofs = col_data.
getIndices().size();
1554 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
1557 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1559 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1561 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
1563 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
1565 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
1567 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
1577 auto v = getVolume();
1578 auto t_w = getFTensor0IntegrationWeight();
1579 auto t_approx_P_adjont_log_du_domega =
1580 getFTensor2FromMat<3, size_symm>(
dataAtPts->adjontPdUdOmegaAtPts);
1582 int row_nb_base_functions = row_data.
getN().size2();
1585 int nb_integration_pts = row_data.
getN().size1();
1587 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1591 for (; rr != row_nb_dofs / 6; ++rr) {
1593 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1594 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1595 double v =
a * t_row_base_fun * t_col_base_fun;
1596 t_m(
L,
k) +=
v * t_approx_P_adjont_log_du_domega(
k,
L);
1603 for (; rr != row_nb_base_functions; ++rr)
1607 ++t_approx_P_adjont_log_du_domega;
1616 int nb_integration_pts = getGaussPts().size2();
1617 int row_nb_dofs = row_data.
getIndices().size();
1618 int col_nb_dofs = col_data.
getIndices().size();
1619 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1622 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1624 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1626 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1636 auto v = getVolume();
1637 auto t_w = getFTensor0IntegrationWeight();
1638 auto t_levi_kirchoff_dP =
1639 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchoffdPAtPts);
1641 int row_nb_base_functions = row_data.
getN().size2();
1644 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1647 for (; rr != row_nb_dofs / 3; ++rr) {
1648 double b =
a * t_row_base_fun;
1650 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1651 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1652 t_m(
k,
i) += b * (t_levi_kirchoff_dP(
i,
l,
k) * t_col_base_fun(
l));
1658 for (; rr != row_nb_base_functions; ++rr) {
1663 ++t_levi_kirchoff_dP;
1671 int nb_integration_pts = getGaussPts().size2();
1672 int row_nb_dofs = row_data.
getIndices().size();
1673 int col_nb_dofs = col_data.
getIndices().size();
1675 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
1677 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
1686 auto v = getVolume();
1687 auto t_w = getFTensor0IntegrationWeight();
1688 auto t_levi_kirchoff_dP =
1689 getFTensor3FromMat<3, 3, 3>(
dataAtPts->leviKirchoffdPAtPts);
1691 int row_nb_base_functions = row_data.
getN().size2();
1694 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1697 for (; rr != row_nb_dofs / 3; ++rr) {
1698 double b =
a * t_row_base_fun;
1699 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1700 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1701 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1702 t_m(
k) += b * (t_levi_kirchoff_dP(
i,
j,
k) * t_col_base_fun(
i,
j));
1709 for (; rr != row_nb_base_functions; ++rr) {
1713 ++t_levi_kirchoff_dP;
1721 int nb_integration_pts = getGaussPts().size2();
1722 int row_nb_dofs = row_data.
getIndices().size();
1723 int col_nb_dofs = col_data.
getIndices().size();
1724 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1727 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1729 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1731 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1742 auto v = getVolume();
1743 auto t_w = getFTensor0IntegrationWeight();
1744 auto t_levi_kirchoff_domega =
1745 getFTensor2FromMat<3, 3>(
dataAtPts->leviKirchoffdOmegaAtPts);
1746 int row_nb_base_functions = row_data.
getN().size2();
1748 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1751 for (; rr != row_nb_dofs / 3; ++rr) {
1752 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1753 const double b =
a * t_row_base_fun;
1755 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1756 t_m(
k,
l) += (b * t_col_base_fun) * t_levi_kirchoff_domega(
k,
l);
1762 for (; rr != row_nb_base_functions; ++rr) {
1766 ++t_levi_kirchoff_domega;
1775 int nb_integration_pts = row_data.
getN().size1();
1776 int row_nb_dofs = row_data.
getIndices().size();
1777 int col_nb_dofs = col_data.
getIndices().size();
1779 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1782 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1784 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1786 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1798 auto v = getVolume();
1799 auto t_w = getFTensor0IntegrationWeight();
1800 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
1801 int row_nb_base_functions = row_data.
getN().size2() / 3;
1804 const double ts_a = getTSa();
1806 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1812 for (; rr != row_nb_dofs / 3; ++rr) {
1815 t_PRT(
i,
k) = t_row_base_fun(
j) * t_h_domega(
i,
j,
k);
1818 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1819 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1820 t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
1828 for (; rr != row_nb_base_functions; ++rr)
1841 int nb_integration_pts = row_data.
getN().size1();
1842 int row_nb_dofs = row_data.
getIndices().size();
1843 int col_nb_dofs = col_data.
getIndices().size();
1845 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1847 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
1857 auto v = getVolume();
1858 auto t_w = getFTensor0IntegrationWeight();
1859 auto t_h_domega = getFTensor3FromMat<3, 3, 3>(
dataAtPts->hdOmegaAtPts);
1860 int row_nb_base_functions = row_data.
getN().size2() / 9;
1863 const double ts_a = getTSa();
1865 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1869 for (; rr != row_nb_dofs; ++rr) {
1872 t_PRT(
k) = t_row_base_fun(
i,
j) * t_h_domega(
i,
j,
k);
1875 auto t_m = get_ftensor2(
K, rr, 0);
1876 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1877 t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
1885 for (; rr != row_nb_base_functions; ++rr)
1898 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1899 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1902 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1907 Tag th_w = create_tag(
"SpatialDisplacement", 3);
1908 Tag th_omega = create_tag(
"Omega", 3);
1909 Tag th_approxP = create_tag(
"Piola1Stress", 9);
1910 Tag th_sigma = create_tag(
"CauchyStress", 9);
1911 Tag th_log_u = create_tag(
"LogSpatialStretch", 9);
1912 Tag th_u = create_tag(
"SpatialStretch", 9);
1913 Tag th_h = create_tag(
"h", 9);
1914 Tag th_X = create_tag(
"X", 3);
1915 Tag th_detF = create_tag(
"detF", 1);
1916 Tag th_angular_momentum = create_tag(
"AngularMomentum", 3);
1918 Tag th_u_eig_vec = create_tag(
"SpatialStretchEigenVec", 9);
1919 Tag th_u_eig_vals = create_tag(
"SpatialStretchEigenVals", 3);
1920 Tag th_traction = create_tag(
"traction", 3);
1922 Tag th_disp = create_tag(
"U", 3);
1923 Tag th_disp_error = create_tag(
"U_ERROR", 1);
1925 auto t_w = getFTensor1FromMat<3>(
dataAtPts->wL2AtPts);
1926 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
1927 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
1929 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStretchTensorAtPts);
1930 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
1931 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1932 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1933 auto t_levi_kirchoff = getFTensor1FromMat<3>(
dataAtPts->leviKirchoffAtPts);
1934 auto t_coords = getFTensor1CoordsAtGaussPts();
1935 auto t_normal = getFTensor1NormalsAtGaussPts();
1936 auto t_disp = getFTensor1FromMat<3>(
dataAtPts->wH1AtPts);
1943 auto set_float_precision = [](
const double x) {
1944 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1951 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1953 v = set_float_precision(
v);
1961 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1964 for (
auto &
a :
v.data())
1965 a = set_float_precision(
a);
1967 &*
v.data().begin());
1975 &
m(0, 0), &
m(0, 1), &
m(0, 2),
1977 &
m(1, 0), &
m(1, 1), &
m(1, 2),
1979 &
m(2, 0), &
m(2, 1), &
m(2, 2));
1981 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1983 t_m(
i,
j) = t_d(
i,
j);
1984 for (
auto &
v :
m.data())
1985 v = set_float_precision(
v);
1987 &*
m.data().begin());
1991 const auto nb_gauss_pts = getGaussPts().size2();
1992 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1995 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
1998 CHKERR save_vec_tag(th_w, t_w, gg);
1999 CHKERR save_vec_tag(th_X, t_coords, gg);
2000 CHKERR save_vec_tag(th_omega, t_omega, gg);
2001 CHKERR save_vec_tag(th_traction, t_traction, gg);
2004 CHKERR save_mat_tag(th_h, t_h, gg);
2007 for (
int ii = 0; ii != 3; ++ii)
2008 for (
int jj = 0; jj != 3; ++jj)
2009 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
2011 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
2014 for (
int ii = 0; ii != 3; ++ii)
2015 for (
int jj = 0; jj != 3; ++jj)
2016 t_u_tmp(ii, jj) = t_u(ii, jj);
2018 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
2019 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
2020 CHKERR save_vec_tag(th_disp, t_disp, gg);
2022 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
2023 CHKERR save_scal_tag(th_disp_error, u_error, gg);
2027 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
2028 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
2032 t_levi(
k) = t_levi_kirchoff(
k);
2036 auto get_eiegn_vector_symmetric = [&](
auto &t_u) {
2039 for (
int ii = 0; ii != 3; ++ii)
2040 for (
int jj = 0; jj != 3; ++jj)
2041 t_m(ii, jj) = t_u(ii, jj);
2044 auto t_eigen_values = getFTensor1FromArray<3>(eigen_values);
2048 &*
m.data().begin());
2050 &*eigen_values.data().begin());
2055 CHKERR get_eiegn_vector_symmetric(t_u);
2077 if (type == MBTET) {
2078 int nb_integration_pts = data.
getN().size1();
2079 auto v = getVolume();
2080 auto t_w = getFTensor0IntegrationWeight();
2081 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
2082 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
2087 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
2088 const double a = t_w *
v;
2090 (*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 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.
VectorBoundedArray< double, 3 > VectorDouble3
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
boost::shared_ptr< BcDispVec > bcDispPtr
MoFEMErrorCode integrate(EntData &data)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< BcRotVec > bcRotPtr
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 integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode integratePiola(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.
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.