14#include <boost/math/constants/constants.hpp>
22double f(
const double v) {
return exp(
v); }
23double d_f(
const double v) {
return exp(
v); }
24double dd_f(
const double v) {
return exp(
v); }
27 typedef typename std::remove_pointer<T>::type
Type;
30 typedef typename std::remove_pointer<T>::type
Type;
45 sqrt(t_omega(
i) * t_omega(
i));
47 constexpr double tol = 1e-18;
48 if (std::abs(angle) <
tol) {
53 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
57 2. * ss_2 * ss_2 / (angle * angle);
58 t_R(
i,
j) +=
a * t_Omega(
i,
j);
59 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
75 sqrt(t_omega(
i) * t_omega(
i));
77 constexpr double tol = 1e-18;
78 if (std::abs(angle) <
tol) {
81 t_diff_R(
i,
j,
k) = FTensor::levi_civita<double>(
i,
l,
k) * t_R(
l,
j);
87 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
94 t_diff_R(
i,
j,
k) = 0;
95 t_diff_R(
i,
j,
k) =
a * FTensor::levi_civita<double>(
i,
j,
k);
97 b * (t_Omega(
i,
l) * FTensor::levi_civita<double>(
l,
j,
k) +
98 FTensor::levi_civita<double>(
i,
l,
k) * t_Omega(
l,
j));
103 (angle * cc - ss) / angle2;
106 (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
107 t_diff_R(
i,
j,
k) += diff_a * t_Omega(
i,
j) * (t_omega(
k) / angle);
109 diff_b * t_Omega(
i,
l) * t_Omega(
l,
j) * (t_omega(
k) / angle);
114template <
typename T1,
typename T2,
typename T3>
118 for (
int ii = 0; ii != 3; ii++)
119 for (
int jj = 0; jj != 3; jj++)
120 eig_vec(ii, jj) = X(ii, jj);
124 int lwork = (3 + 2) * 3;
125 std::array<
double, (3 + 2) * 3> work;
127 if (
lapack_dsyev(
'V',
'U',
n, &(eig_vec(0, 0)), lda, &eig(0), work.data(),
130 "The algorithm failed to compute eigenvalues.");
134inline bool is_eq(
const double &
a,
const double &b) {
135 constexpr double eps = 1e-8;
136 return abs(
a - b) <
eps;
140 return distance(&ec(0), unique(&ec(0), &ec(0) + 3,
is_eq));
143template <
typename T1,
typename T2>
147 if (
is_eq(eig(0), eig(1))) {
149 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
150 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
151 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
154 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
157 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
158 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
159 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
162 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
172 int nb_integration_pts = data.
getN().size1();
178 dataAtPts->hAtPts.resize(9, nb_integration_pts,
false);
179 dataAtPts->rotMatAtPts.resize(9, nb_integration_pts,
false);
180 dataAtPts->diffRotMatAtPts.resize(27, nb_integration_pts,
false);
181 dataAtPts->streachTensorAtPts.resize(6, nb_integration_pts,
false);
182 dataAtPts->diffStreachTensorAtPts.resize(36, nb_integration_pts,
false);
183 dataAtPts->eigenVals.resize(3, nb_integration_pts,
false);
184 dataAtPts->eigenVecs.resize(9, nb_integration_pts,
false);
185 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
187 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
188 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
189 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
192 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachTensorAtPts);
193 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
196 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
197 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
198 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
201 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
206 t_diff_R(
i,
j,
k) = t0_diff(
i,
j,
k);
207 t_R(
i,
j) = t0(
i,
j);
213 t_eigen_vals(
i) = eig(
i);
214 t_eigen_vecs(
i,
j) = eigen_vec(
i,
j);
224 for (
int ii = 0; ii != 3; ++ii) {
225 for (
int jj = ii; jj != 3; ++jj) {
226 for (
int kk = 0; kk != 3; ++kk) {
227 for (
int ll = 0; ll < kk; ++ll) {
228 t_diff_u_tmp(ii, jj, kk, ll) *= 2;
234 t_u(
i,
j) = t_u_tmp(
i,
j);
235 t_diff_u(
i,
j,
k,
l) = t_diff_u_tmp(
i,
j,
k,
l);
236 t_h(
i,
j) = t_R(
i,
k) * t_u(
k,
j);
255 int nb_integration_pts = data.
getN().size1();
256 auto v = getVolume();
257 auto t_w = getFTensor0IntegrationWeight();
258 auto t_div_P = getFTensor1FromMat<3>(
dataAtPts->divPAtPts);
259 auto t_s_dot_w = getFTensor1FromMat<3>(
dataAtPts->wDotAtPts);
260 if (
dataAtPts->wDotDotAtPts.size1() == 0 &&
261 dataAtPts->wDotDotAtPts.size2() != nb_integration_pts) {
262 dataAtPts->wDotDotAtPts.resize(3, nb_integration_pts);
265 auto t_s_dot_dot_w = getFTensor1FromMat<3>(
dataAtPts->wDotDotAtPts);
267 int nb_base_functions = data.
getN().size2();
271 auto get_ftensor1 = [](
auto &
v) {
276 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
278 auto t_nf = get_ftensor1(
nF);
280 for (; bb != nb_dofs / 3; ++bb) {
281 t_nf(
i) +=
a * t_row_base_fun * t_div_P(
i);
282 t_nf(
i) -=
a * t_row_base_fun *
alphaW * t_s_dot_w(
i);
283 t_nf(
i) -=
a * t_row_base_fun *
alphaRho * t_s_dot_dot_w(
i);
287 for (; bb != nb_base_functions; ++bb)
301 int nb_integration_pts = data.
getN().size1();
302 auto v = getVolume();
303 auto t_w = getFTensor0IntegrationWeight();
304 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
305 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
307 int nb_base_functions = data.
getN().size2();
314 auto get_ftensor1 = [](
auto &
v) {
319 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
321 auto t_nf = get_ftensor1(
nF);
323 t_PRT(
i,
m) = t_approx_P(
i,
j) * t_R(
m,
j);
325 t_leviPRT(
k) = levi_civita(
i,
m,
k) * t_PRT(
i,
m);
327 for (; bb != nb_dofs / 3; ++bb) {
328 t_nf(
k) +=
a * t_row_base_fun * t_leviPRT(
k);
332 for (; bb != nb_base_functions; ++bb)
344 int nb_integration_pts = data.
getN().size1();
345 auto v = getVolume();
346 auto t_w = getFTensor0IntegrationWeight();
347 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
348 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
349 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
351 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
353 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
359 auto get_ftensor2 = [](
auto &
v) {
361 &
v[0], &
v[1], &
v[2], &
v[3], &
v[4], &
v[5]);
364 int nb_base_functions = data.
getN().size2();
366 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
368 auto t_nf = get_ftensor2(
nF);
371 t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
373 t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
376 t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
380 a * (t_diff_u(
k,
l,
i,
j) * (t_deltaP(
k,
l) - t_viscous_stress(
k,
l)));
383 for (; bb != nb_dofs / 6; ++bb) {
385 t_nf(
i,
j) += t_row_base_fun * t1(
i,
j);
390 for (; bb != nb_base_functions; ++bb)
406 int nb_integration_pts = data.
getN().size1();
407 auto v = getVolume();
408 auto t_w = getFTensor0IntegrationWeight();
409 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
410 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
411 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
413 int nb_base_functions = data.
getN().size2() / 3;
420 auto get_ftensor1 = [](
auto &
v) {
425 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
427 auto t_nf = get_ftensor1(
nF);
433 t_residuum_omega(
i,
j) =
434 (levi_civita(
i,
m,
k) * t_omega_dot(
k)) * t_R(
m,
j);
437 for (; bb != nb_dofs / 3; ++bb) {
438 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum(
i,
j);
439 t_nf(
i) +=
a * t_row_base_fun(
j) * t_residuum_omega(
i,
j);
445 for (; bb != nb_base_functions; ++bb)
460 int nb_integration_pts = data.
getN().size1();
461 auto v = getVolume();
462 auto t_w = getFTensor0IntegrationWeight();
463 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
464 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
465 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
467 int nb_base_functions = data.
getN().size2() / 9;
474 auto get_ftensor0 = [](
auto &
v) {
478 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
480 auto t_nf = get_ftensor0(
nF);
486 t_residuum_omega(
i,
j) =
487 (levi_civita(
i,
m,
k) * t_omega_dot(
k)) * t_R(
m,
j);
490 for (; bb != nb_dofs; ++bb) {
491 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum(
i,
j);
492 t_nf +=
a * t_row_base_fun(
i,
j) * t_residuum_omega(
i,
j);
496 for (; bb != nb_base_functions; ++bb) {
511 int nb_integration_pts = data.
getN().size1();
512 auto v = getVolume();
513 auto t_w = getFTensor0IntegrationWeight();
514 auto t_s_w = getFTensor1FromMat<3>(
dataAtPts->wAtPts);
515 int nb_base_functions = data.
getN().size2() / 3;
518 auto get_ftensor1 = [](
auto &
v) {
523 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
525 auto t_nf = get_ftensor1(
nF);
527 for (; bb != nb_dofs / 3; ++bb) {
528 double div_row_base = t_row_diff_base_fun(
i,
i);
529 t_nf(
i) +=
a * div_row_base * t_s_w(
i);
531 ++t_row_diff_base_fun;
533 for (; bb != nb_base_functions; ++bb) {
534 ++t_row_diff_base_fun;
550 if (bc.faces.find(fe_ent) != bc.faces.end()) {
552 int nb_integration_pts = data.
getN().size1();
553 auto t_normal = getFTensor1Normal();
554 auto t_w = getFTensor0IntegrationWeight();
555 int nb_base_functions = data.
getN().size2() / 3;
560 auto get_ftensor1 = [](
auto &
v) {
567 t_bc_disp(
i) *= getFEMethod()->ts_t;
569 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
571 t_bc_res(
i) = t_bc_disp(
i);
572 auto t_nf = get_ftensor1(
nF);
574 for (; bb != nb_dofs / 3; ++bb) {
576 t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_bc_res(
i) * 0.5;
580 for (; bb != nb_base_functions; ++bb)
597 if (bc.faces.find(fe_ent) != bc.faces.end()) {
599 int nb_integration_pts = data.
getN().size1();
600 auto t_normal = getFTensor1Normal();
601 auto t_w = getFTensor0IntegrationWeight();
603 int nb_base_functions = data.
getN().size2() / 3;
608 auto get_ftensor1 = [](
auto &
v) {
616 double theta = bc.theta;
617 theta *= getFEMethod()->ts_t;
620 const double a = sqrt(t_normal(
i) * t_normal(
i));
621 t_omega(
i) = t_normal(
i) * (theta /
a);
623 auto t_coords = getFTensor1CoordsAtGaussPts();
625 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
627 t_delta(
i) = t_center(
i) - t_coords(
i);
629 t_disp(
i) = t_delta(
i) - t_R(
i,
j) * t_delta(
j);
631 auto t_nf = get_ftensor1(
nF);
633 for (; bb != nb_dofs / 3; ++bb) {
634 t_nf(
i) -= t_w * (t_row_base_fun(
j) * t_normal(
j)) * t_disp(
i) * 0.5;
638 for (; bb != nb_base_functions; ++bb)
653 int operator()(
int p_row,
int p_col,
int p_data)
const {
654 return 2 * (p_data + 1);
658 if (
ts_ctx == CTX_TSSETIFUNCTION) {
670 this->getCacheWeakPtr());
682 auto t_normal = getFTensor1Normal();
683 const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
685 t_unit_normal(
i) = t_normal(
i) / nrm2;
687 int nb_integration_pts = data.
getN().size1();
688 int nb_base_functions = data.
getN().size2() / 3;
689 double ts_t = getFEMethod()->ts_t;
691 auto integrate_matrix = [&]() {
694 auto t_w = getFTensor0IntegrationWeight();
695 matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
699 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
702 for (; rr != nb_dofs / 3; ++rr) {
703 const double a = t_row_base_fun(
i) * t_unit_normal(
i);
705 for (
int cc = 0; cc != nb_dofs / 3; ++cc) {
706 const double b = t_col_base_fun(
i) * t_unit_normal(
i);
707 matPP(rr, cc) += t_w *
a * b;
713 for (; rr != nb_base_functions; ++rr)
722 auto integrate_rhs = [&](
auto &bc) {
725 auto t_w = getFTensor0IntegrationWeight();
726 vecPv.resize(3, nb_dofs / 3,
false);
730 double ts_t = getFEMethod()->ts_t;
732 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
734 for (; rr != nb_dofs / 3; ++rr) {
735 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
736 for (
int dd = 0; dd != 3; ++dd)
738 vecPv(dd, rr) +=
t * bc.vals[dd];
741 for (; rr != nb_base_functions; ++rr)
748 auto integrate_rhs_cook = [&](
auto &bc) {
751 vecPv.resize(3, nb_dofs / 3,
false);
754 auto t_w = getFTensor0IntegrationWeight();
756 auto t_coords = getFTensor1CoordsAtGaussPts();
758 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
760 auto calc_tau = [](
double y) {
763 return -y * (y - 1) / 0.25;
766 const double tau = calc_tau(t_coords(1));
769 for (; rr != nb_dofs / 3; ++rr) {
770 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
771 for (
int dd = 0; dd != 3; ++dd)
773 vecPv(dd, rr) +=
t * tau * bc.vals[dd];
777 for (; rr != nb_base_functions; ++rr)
788 if (bc.faces.find(fe_ent) != bc.faces.end()) {
793 CHKERR integrate_matrix();
794 if (bc.blockName.compare(20, 4,
"COOK") == 0)
795 CHKERR integrate_rhs_cook(bc);
801 nb_dofs / 3, &*
vecPv.data().begin(), nb_dofs / 3);
804 "The leading minor of order %i is "
805 "not positive; definite;\nthe "
806 "solution could not be computed",
809 for (
int dd = 0; dd != 3; ++dd)
811 for (
int rr = 0; rr != nb_dofs / 3; ++rr)
823 int nb_integration_pts = row_data.
getN().size1();
824 int row_nb_dofs = row_data.
getIndices().size();
825 int col_nb_dofs = col_data.
getIndices().size();
826 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
828 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2));
831 auto v = getVolume();
832 auto t_w = getFTensor0IntegrationWeight();
833 int row_nb_base_functions = row_data.
getN().size2();
835 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
838 for (; rr != row_nb_dofs / 3; ++rr) {
840 auto t_m = get_ftensor1(
K, 3 * rr, 0);
841 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
842 double div_col_base = t_col_diff_base_fun(
i,
i);
843 t_m(
i) +=
a * t_row_base_fun * div_col_base;
845 ++t_col_diff_base_fun;
849 for (; rr != row_nb_base_functions; ++rr)
860 if (
alphaW < std::numeric_limits<double>::epsilon() &&
861 alphaRho < std::numeric_limits<double>::epsilon())
864 const int nb_integration_pts = row_data.
getN().size1();
865 const int row_nb_dofs = row_data.
getIndices().size();
866 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
868 &
m(r + 0,
c + 0), &
m(r + 1,
c + 1), &
m(r + 2,
c + 2)
874 auto v = getVolume();
875 auto t_w = getFTensor0IntegrationWeight();
877 int row_nb_base_functions = row_data.
getN().size2();
880 double ts_scale =
alphaW * getTSa();
881 if (std::abs(
alphaRho) > std::numeric_limits<double>::epsilon())
884 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
885 double a =
v * t_w * ts_scale;
888 for (; rr != row_nb_dofs / 3; ++rr) {
891 auto t_m = get_ftensor1(
K, 3 * rr, 0);
892 for (
int cc = 0; cc != row_nb_dofs / 3; ++cc) {
893 const double b =
a * t_row_base_fun * t_col_base_fun;
902 for (; rr != row_nb_base_functions; ++rr)
914 int nb_integration_pts = row_data.
getN().size1();
915 int row_nb_dofs = row_data.
getIndices().size();
916 int col_nb_dofs = col_data.
getIndices().size();
917 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
920 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
922 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
924 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
926 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
928 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
930 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2));
940 auto v = getVolume();
941 auto t_w = getFTensor0IntegrationWeight();
942 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
945 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
947 int row_nb_base_functions = row_data.
getN().size2();
949 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
953 for (; rr != row_nb_dofs / 6; ++rr) {
956 auto t_m = get_ftensor3(
K, 6 * rr, 0);
958 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
961 t_RTP_dP(
i,
j,
k) = t_R(
k,
i) * t_col_base_fun(
j) * t_row_base_fun;
963 for (
int kk = 0; kk != 3; ++kk)
964 for (
int ii = 0; ii != 3; ++ii)
965 for (
int jj = ii; jj != 3; ++jj) {
966 for (
int mm = 0; mm != 3; ++mm)
967 for (
int nn = 0; nn != 3; ++nn)
969 a * t_diff_u(mm, nn, ii, jj) * t_RTP_dP(mm, nn, kk);
978 for (; rr != row_nb_base_functions; ++rr)
990 int nb_integration_pts = row_data.
getN().size1();
991 int row_nb_dofs = row_data.
getIndices().size();
992 int col_nb_dofs = col_data.
getIndices().size();
993 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
995 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c), &
m(r + 3,
c), &
m(r + 4,
c),
1005 auto v = getVolume();
1006 auto t_w = getFTensor0IntegrationWeight();
1007 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1009 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
1011 int row_nb_base_functions = row_data.
getN().size2();
1013 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1017 for (; rr != row_nb_dofs / 6; ++rr) {
1019 auto t_m = get_ftensor2(
K, 6 * rr, 0);
1021 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1022 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1025 t_RTP_dBubble(
i,
j) =
1026 a * t_row_base_fun * t_R(
k,
i) * t_col_base_fun(
k,
j);
1027 t_m(
i,
j) += t_diff_u(
m,
n,
i,
j) * t_RTP_dBubble(
m,
n);
1035 for (; rr != row_nb_base_functions; ++rr)
1048 int nb_integration_pts = row_data.
getN().size1();
1049 int row_nb_dofs = row_data.
getIndices().size();
1050 int col_nb_dofs = col_data.
getIndices().size();
1051 auto get_ftensor4 = [](
MatrixDouble &
m,
const int r,
const int c) {
1054 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1055 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1057 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1058 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1060 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1061 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1063 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1064 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1066 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1067 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1069 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1070 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1085 auto v = getVolume();
1086 auto t_w = getFTensor0IntegrationWeight();
1090 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1091 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1092 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
1094 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
1095 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1097 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
1098 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1099 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1107 t_one4_symmetric(
i,
j,
k,
l) = 0;
1108 for (
auto ii : {0, 1, 2})
1109 for (
auto jj : {0, 1, 2})
1110 for (
auto kk : {0, 1, 2})
1111 for (
auto ll : {0, 1, 2}) {
1114 t_one4_symmetric(ii, jj, kk, ll) =
1115 t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1117 t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1120 int row_nb_base_functions = row_data.
getN().size2();
1125 for (
auto ii : {0, 1, 2})
1128 const double ts_a = getTSa();
1129 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1133 t_P_dh(
i,
j,
N0,
k) = t_P_dh0(
i,
j,
k);
1134 t_P_dh(
i,
j,
N1,
k) = t_P_dh1(
i,
j,
k);
1135 t_P_dh(
i,
j,
N2,
k) = t_P_dh2(
i,
j,
k);
1138 t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
1140 t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
1143 t_dot_u(
i,
j) = t_u(
i,
k) * t_dot_log_u(
k,
j);
1146 t_stress_diff(
i,
j,
k,
l) = 0;
1147 for (
int ii = 0; ii != 3; ++ii)
1148 for (
int jj = 0; jj != 3; ++jj)
1149 for (
int kk = 0; kk != 3; ++kk)
1150 for (
int ll = 0; ll != 3; ++ll)
1151 for (
int mm = 0; mm != 3; ++mm)
1152 for (
int nn = 0; nn != 3; ++nn)
1153 t_stress_diff(ii, jj, mm, nn) -=
1154 t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1157 t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
1159 t_viscous_stress_diff(
i,
j,
k,
l) = t_one4_symmetric(
i,
j,
k,
l);
1160 t_viscous_stress_diff(
i,
j,
k,
l) *= (ts_a *
alphaU);
1163 t_deltaP2(
i,
j) = t_deltaP(
i,
j) - t_viscous_stress(
i,
j);
1166 t_eigen_vals, t_eigen_vecs,
f,
d_f,
dd_f, t_deltaP2, nbUniq[gg]);
1167 for (
int ii = 0; ii != 3; ++ii) {
1168 for (
int jj = 0; jj < ii; ++jj) {
1169 for (
int kk = 0; kk != 3; ++kk) {
1170 for (
int ll = 0; ll != kk; ++ll) {
1171 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1176 for (
int ii = 0; ii != 3; ++ii) {
1177 for (
int jj = ii; jj != 3; ++jj) {
1178 for (
int kk = 0; kk != 3; ++kk) {
1179 for (
int ll = 0; ll < kk; ++ll) {
1180 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1187 for (; rr != row_nb_dofs / 6; ++rr) {
1190 auto t_m = get_ftensor4(
K, 6 * rr, 0);
1192 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1194 const double b =
a * t_row_base_fun * t_col_base_fun;
1196 for (
int ii = 0; ii != 3; ++ii)
1197 for (
int jj = ii; jj != 3; ++jj)
1198 for (
int kk = 0; kk != 3; ++kk)
1199 for (
int ll = kk; ll != 3; ++ll)
1200 for (
int mm = 0; mm != 3; ++mm)
1201 for (
int nn = 0; nn != 3; ++nn)
1202 t_m(ii, jj, kk, ll) +=
1203 b * t_diff_u(mm, nn, ii, jj) *
1204 (t_stress_diff(mm, nn, kk, ll) -
1205 t_viscous_stress_diff(mm, nn, kk, ll));
1207 t_m(
i,
j,
k,
l) += b * t_diff2_uP2(
i,
j,
k,
l);
1215 for (; rr != row_nb_base_functions; ++rr)
1237 int row_nb_dofs = row_data.
getIndices().size();
1238 int col_nb_dofs = col_data.
getIndices().size();
1239 auto get_ftensor3 = [](
MatrixDouble &
m,
const int r,
const int c) {
1242 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1244 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1246 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2),
1248 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2),
1250 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2),
1252 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2)
1261 auto v = getVolume();
1262 auto t_w = getFTensor0IntegrationWeight();
1263 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1266 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
1268 int row_nb_base_functions = row_data.
getN().size2();
1271 int nb_integration_pts = row_data.
getN().size1();
1273 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1277 t_RT_domega_P(
i,
j,
m) = t_diff_R(
k,
i,
m) * t_P(
k,
j);
1280 for (
int ii = 0; ii != 3; ++ii)
1281 for (
int jj = ii; jj != 3; ++jj)
1282 for (
int mm = 0; mm != 3; ++mm)
1283 for (
int nn = 0; nn != 3; ++nn)
1284 for (
int kk = 0; kk != 3; ++kk)
1286 t_diff_u(mm, nn, ii, jj) * t_RT_domega_P(mm, nn, kk);
1289 for (; rr != row_nb_dofs / 6; ++rr) {
1291 auto t_m = get_ftensor3(
K, 6 * rr, 0);
1292 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1294 double v =
a * t_row_base_fun * t_col_base_fun;
1304 for (; rr != row_nb_base_functions; ++rr)
1319 int nb_integration_pts = row_data.
getN().size1();
1320 int row_nb_dofs = row_data.
getIndices().size();
1321 int col_nb_dofs = col_data.
getIndices().size();
1322 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1325 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1327 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1329 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1340 auto v = getVolume();
1341 auto t_w = getFTensor0IntegrationWeight();
1342 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1344 int row_nb_base_functions = row_data.
getN().size2();
1347 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1351 for (; rr != row_nb_dofs / 3; ++rr) {
1353 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1356 (
a * t_row_base_fun) * (levi_civita(
i,
m,
k) * t_R(
m,
j));
1357 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1359 t_m(
k,
i) += t_levi_R(
i,
j,
k) * t_col_base_fun(
j);
1367 for (; rr != row_nb_base_functions; ++rr)
1378 int nb_integration_pts = row_data.
getN().size1();
1379 int row_nb_dofs = row_data.
getIndices().size();
1380 int col_nb_dofs = col_data.
getIndices().size();
1381 auto get_ftensor1 = [](
MatrixDouble &
m,
const int r,
const int c) {
1383 &
m(r + 0,
c), &
m(r + 1,
c), &
m(r + 2,
c));
1391 auto v = getVolume();
1392 auto t_w = getFTensor0IntegrationWeight();
1393 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1395 int row_nb_base_functions = row_data.
getN().size2();
1398 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1402 for (; rr != row_nb_dofs / 3; ++rr) {
1404 auto t_col_base_fun = col_data.
getFTensor2N<3, 3>(gg, 0);
1405 auto t_m = get_ftensor1(
K, 3 * rr, 0);
1409 (
a * t_row_base_fun) * (levi_civita(
i,
m,
k) * t_R(
m,
j));
1411 for (
int cc = 0; cc != col_nb_dofs; ++cc) {
1413 t_m(
k) += t_levi_r(
i,
j,
k) * t_col_base_fun(
i,
j);
1421 for (; rr != row_nb_base_functions; ++rr)
1432 int nb_integration_pts = row_data.
getN().size1();
1433 int row_nb_dofs = row_data.
getIndices().size();
1434 int col_nb_dofs = col_data.
getIndices().size();
1435 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),
1440 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1442 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1453 auto v = getVolume();
1454 auto t_w = getFTensor0IntegrationWeight();
1455 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1456 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1459 int row_nb_base_functions = row_data.
getN().size2();
1462 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1466 t_PRT(
i,
m) = t_approx_P(
i,
j) * t_R(
m,
j);
1469 t_leviPRT_domega(
n,
l) =
1470 levi_civita(
i,
j,
n) * (t_approx_P(
i,
k) * t_diff_R(
j,
k,
l));
1473 for (; rr != row_nb_dofs / 3; ++rr) {
1476 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1477 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1479 const double b =
a * t_row_base_fun * t_col_base_fun;
1481 t_m(
k,
m) += b * t_leviPRT_domega(
k,
m);
1490 for (; rr != row_nb_base_functions; ++rr)
1504 int nb_integration_pts = row_data.
getN().size1();
1505 int row_nb_dofs = row_data.
getIndices().size();
1506 int col_nb_dofs = col_data.
getIndices().size();
1508 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1511 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2),
1513 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2),
1515 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2)
1527 auto v = getVolume();
1528 auto t_w = getFTensor0IntegrationWeight();
1529 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1531 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1532 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
1533 int row_nb_base_functions = row_data.
getN().size2() / 3;
1536 const double ts_a = getTSa();
1538 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1544 for (; rr != row_nb_dofs / 3; ++rr) {
1547 t_PRT(
i,
k) = t_row_base_fun(
j) * (t_diff_R(
i,
l,
k) * t_u(
l,
j));
1550 t_levi1(
i,
k) = (levi_civita(
i,
m,
n) * t_omega_dot(
n)) *
1551 (t_row_base_fun(
j) * t_diff_R(
m,
j,
k));
1552 t_levi2(
i,
k) = levi_civita(
i,
m,
k) * (t_row_base_fun(
j) * t_R(
m,
j));
1555 auto t_m = get_ftensor2(
K, 3 * rr, 0);
1556 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1557 t_m(
i,
j) += (
a * t_col_base_fun) * t_PRT(
i,
j);
1558 t_m(
i,
j) += (
a * t_col_base_fun) * t_levi1(
i,
j);
1559 t_m(
i,
j) += ((
a * ts_a) * t_col_base_fun) * t_levi2(
i,
j);
1567 for (; rr != row_nb_base_functions; ++rr)
1583 int nb_integration_pts = row_data.
getN().size1();
1584 int row_nb_dofs = row_data.
getIndices().size();
1585 int col_nb_dofs = col_data.
getIndices().size();
1587 auto get_ftensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
1589 &
m(r,
c + 0), &
m(r,
c + 1), &
m(r,
c + 2));
1599 auto v = getVolume();
1600 auto t_w = getFTensor0IntegrationWeight();
1601 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1603 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1604 auto t_omega_dot = getFTensor1FromMat<3>(
dataAtPts->rotAxisDotAtPts);
1605 int row_nb_base_functions = row_data.
getN().size2() / 9;
1608 const double ts_a = getTSa();
1610 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1614 for (; rr != row_nb_dofs; ++rr) {
1617 t_PRT(
k) = t_row_base_fun(
i,
j) * (t_diff_R(
i,
l,
k) * t_u(
l,
j));
1620 t_levi_omega(
i,
m) = levi_civita(
i,
m,
n) * t_omega_dot(
n);
1622 t_levi_omegaDiffR(
i,
j,
k) = t_levi_omega(
i,
m) * t_diff_R(
m,
j,
k);
1624 t_levi_omegaR(
i,
j,
k) = levi_civita(
i,
m,
k) * t_R(
m,
j);
1627 t_levi1(
k) = t_row_base_fun(
i,
j) * t_levi_omegaDiffR(
i,
j,
k);
1628 t_levi2(
k) = t_row_base_fun(
i,
j) * t_levi_omegaR(
i,
j,
k);
1631 auto t_m = get_ftensor2(
K, rr, 0);
1632 for (
int cc = 0; cc != col_nb_dofs / 3; ++cc) {
1633 t_m(
j) += (
a * t_col_base_fun) * t_PRT(
j);
1634 t_m(
j) += (
a * t_col_base_fun) * t_levi1(
j);
1635 t_m(
j) += ((
a * ts_a) * t_col_base_fun) * t_levi2(
j);
1643 for (; rr != row_nb_base_functions; ++rr)
1658 if (type != MBVERTEX)
1661 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1662 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1665 th, MB_TAG_CREAT | MB_TAG_SPARSE,
1670 Tag th_w = create_tag(
"SpatialDisplacement", 3);
1671 Tag th_omega = create_tag(
"Omega", 3);
1672 Tag th_approxP = create_tag(
"Piola1Stress", 9);
1673 Tag th_sigma = create_tag(
"CauchyStress", 9);
1674 Tag th_log_u = create_tag(
"LogSpatialStreach", 9);
1675 Tag th_u = create_tag(
"SpatialStreach", 9);
1676 Tag th_h = create_tag(
"h", 9);
1677 Tag th_X = create_tag(
"X", 3);
1678 Tag th_detF = create_tag(
"detF", 1);
1679 Tag th_angular_momentum = create_tag(
"AngularMomentum", 3);
1681 Tag th_u_eig_vec = create_tag(
"SpatialStreachEigenVec", 9);
1682 Tag th_u_eig_vals = create_tag(
"SpatialStreachEigenVals", 3);
1684 int nb_gauss_pts = data.
getN().size1();
1685 if (
mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
1687 "Nb. of integration points is not equal to number points on "
1688 "post-processing mesh");
1691 auto t_w = getFTensor1FromMat<3>(
dataAtPts->wAtPts);
1692 auto t_omega = getFTensor1FromMat<3>(
dataAtPts->rotAxisAtPts);
1693 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
1695 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachTensorAtPts);
1696 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1697 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1698 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1699 auto t_coords = getFTensor1CoordsAtGaussPts();
1709 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1713 &*
v.data().begin());
1719 &
m(0, 0), &
m(0, 1), &
m(0, 2),
1721 &
m(1, 0), &
m(1, 1), &
m(1, 2),
1723 &
m(2, 0), &
m(2, 1), &
m(2, 2));
1725 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1727 t_m(
i,
j) = t_d(
i,
j);
1729 &*
m.data().begin());
1733 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
1736 CHKERR save_vec_tag(th_w, t_w, gg);
1737 CHKERR save_vec_tag(th_X, t_coords, gg);
1738 CHKERR save_vec_tag(th_omega, t_omega, gg);
1741 CHKERR save_mat_tag(th_h, t_h, gg);
1744 for (
int ii = 0; ii != 3; ++ii)
1745 for (
int jj = 0; jj != 3; ++jj)
1746 t_log_u_tmp(ii, jj) = t_log_u(ii, jj);
1748 CHKERR save_mat_tag(th_log_u, t_log_u_tmp, gg);
1751 for (
int ii = 0; ii != 3; ++ii)
1752 for (
int jj = 0; jj != 3; ++jj)
1753 t_u_tmp(ii, jj) = t_u(ii, jj);
1755 CHKERR save_mat_tag(th_u, t_u_tmp, gg);
1756 CHKERR save_mat_tag(th_approxP, t_approx_P, gg);
1760 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
1761 CHKERR save_mat_tag(th_sigma, t_cauchy, gg);
1765 t_PhT(
i,
k) = t_approx_P(
i,
j) * t_R(
k,
j);
1767 t_leviPRT(
k) = levi_civita(
i,
l,
k) * t_PhT(
i,
l);
1772 auto get_eiegn_vector_symmetric = [&](
auto &t_u) {
1775 for (
int ii = 0; ii != 3; ++ii)
1776 for (
int jj = 0; jj != 3; ++jj)
1777 t_m(ii, jj) = t_u(ii, jj);
1783 int n = 3, lda = 3, info, lwork = -1;
1786 &(eigen_values.data()[0]), &wkopt, lwork);
1789 "is something wrong with lapack_dsyev info = %d", info);
1793 &(eigen_values.data()[0]), work, lwork);
1796 "is something wrong with lapack_dsyev info = %d", info);
1799 &*
m.data().begin());
1801 &*eigen_values.data().begin());
1806 CHKERR get_eiegn_vector_symmetric(t_u);
1825 if (type == MBTET) {
1826 int nb_integration_pts = data.
getN().size1();
1827 auto v = getVolume();
1828 auto t_w = getFTensor0IntegrationWeight();
1829 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1830 auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
1835 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1836 const double a = t_w *
v;
1837 (*energy) +=
a * t_P(
i,
J) * t_h(
i,
J);
Eshelbian plasticity interface.
#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
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
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)
static FTensor::Tensor2< typename TensorTypeExtractor< T >::Type, 3, 3 > get_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
MoFEMErrorCode get_eigen_val_and_proj_lapack(T1 &X, T2 &eig, T3 &eig_vec)
bool is_eq(const double &a, const double &b)
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
void sort_eigen_vals(T1 &eig, T2 &eigen_vec)
double dd_f(const double v)
static FTensor::Tensor3< typename TensorTypeExtractor< T >::Type, 3, 3, 3 > get_diff_rotation_form_vector(FTensor::Tensor1< T, 3 > &t_omega)
double d_f(const double v)
Tensors class implemented by Walter Landry.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
implementation of Data Operators for Forces and Sources
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
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 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 integrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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.