33 boost::shared_ptr<MatrixDouble>
mtD;
43 boost::shared_ptr<MatrixDouble>
matC;
44 boost::shared_ptr<VectorDouble>
detF;
91 Tensor2<double, 3, 3> non_symm;
104template <
typename T>
inline auto to_symm(
T &non_symm) {
105 Tensor2_symmetric<double, 3> symm;
118inline bool is_eq(
const double &a,
const double &b) {
119 constexpr double eps = 1e-12;
120 return abs(
a - b) <
eps;
122inline bool is_diff(
const double &a,
const double &b) {
123 constexpr double eps = 1e-12;
124 return abs(
a - b) >
eps;
129 Ddg<double, 3, 3> loc_tens;
176 Tensor4<double, 3, 3, 3, 3> tens4;
177 for (
int ii = 0; ii != 3; ++ii)
178 for (
int jj = 0; jj != 3; ++jj)
179 for (
int kk = 0; kk != 3; ++kk)
180 for (
int ll = 0; ll != 3; ++ll)
181 tens4(ii, jj, kk, ll) = ddg(ii, jj, kk, ll);
185template <
typename T>
inline bool is_ddg(Tensor4<T, 3, 3, 3, 3> tensor) {
187 for (
int ii = 0; ii != 3; ++ii)
188 for (
int jj = 0; jj != 3; ++jj)
189 for (
int kk = 0; kk != 3; ++kk)
190 for (
int ll : {0, 1, 2}) {
192 abs(test_ddg(ii, jj, kk, ll) - tensor(ii, jj, kk, ll));
200template <
typename T>
inline auto exp_map(
const Tensor2_symmetric<T, 3> &X) {
201 Tensor2<double, 3, 3> exp_tensor(1., 0., 0., 0., 1., 0., 0., 0., 1.);
202 Tensor2<double, 3, 3> X_nt;
204 constexpr int max_it = 1000;
205 constexpr double eps = 1e-12;
211 exp_tensor(
i,
j) += Xcopy(
i,
j);
213 for (
int n = 2;
n != max_it; ++
n) {
215 X_nt(
i,
j) = X_n(
i,
j);
216 X_n(
i,
k) = X_nt(
i,
j) * Xcopy(
j,
k);
217 const double inv_nf = 1. / fac_n;
218 exp_tensor(
i,
j) += inv_nf * X_n(
i,
j);
219 const double e = inv_nf * std::abs(sqrt(X_n(
i,
j) * X_n(
i,
j)));
236 const double angle = sqrt(t_omega(
i) * t_omega(
i));
237 if (std::abs(angle) < 1e-18)
241 t_Omega(
i,
j) = FTensor::levi_civita<double>(
i,
j,
k) * t_omega(
k);
242 const double a = sin(angle) / angle;
243 const double ss_2 = sin(angle / 2.);
244 const double b = 2. * ss_2 * ss_2 / (angle * angle);
245 t_R(
i,
j) +=
a * t_Omega(
i,
j);
246 t_R(
i,
j) += b * t_Omega(
i,
k) * t_Omega(
k,
j);
252 my_c_omega(
i) = t1_omega(
i) * tt;
253 auto t_R = get_rotation(my_c_omega);
255 return get_rotation(my_c_omega);
261 constexpr int max_it = 1000;
262 constexpr double eps = 1e-12;
265 log_tensor(
i,
j) -= kronecker_delta(
i,
j);
269 Tensor2<double, 3, 3> X_nt_tmp;
271 for (
int n = 2;
n != max_it; ++
n) {
273 X_nt_tmp(
i,
j) = X_n(
i,
j);
274 X_n(
i,
k) = X_nt_tmp(
i,
j) * Xcopy(
j,
k);
275 const double inv_n = 1. /
n;
276 log_tensor(
i,
j) += inv_n * pow(-1,
n + 1) * X_n(
i,
j);
277 const double e = inv_n * std::abs(sqrt(X_n(
i,
j) * X_n(
i,
j)));
286 Tensor2<double, 3, 3> X_nt_tmp;
287 Tensor4<double, 3, 3, 3, 3> diff_log_map;
288 Tensor4<double, 3, 3, 3, 3> diff_log_map_tmp;
290 constexpr int max_it = 1000;
291 constexpr double eps = 1e-12;
293 auto get_Xn = [&](
auto &log_tensor) {
294 vector<Tensor2<double, 3, 3>> Xn;
296 auto &cu_Xn = Xn.back();
297 cu_Xn(
i,
j) = kronecker_delta(
i,
j);
298 for (
int n = 2;
n != max_it; ++
n) {
299 const double cof = 1. /
n;
300 auto X_nt_tmp = Xn.back();
302 auto &c_Xn = Xn.back();
303 c_Xn(
i,
k) = X_nt_tmp(
i,
j) * log_tensor(
j,
k);
304 const double e = cof * std::abs(sqrt(c_Xn(
i,
j) * c_Xn(
i,
j)));
312 log_tensor(
i,
j) -= kronecker_delta(
i,
j);
313 diff_log_map(
i,
j,
k,
l) = 0.;
315 auto Xn_vec = get_Xn(log_tensor);
317 for (
int n = 1;
n != Xn_vec.size() + 1; ++
n) {
318 const double cof = pow(-1,
n + 1) /
n;
319 for (
int m = 1;
m !=
n + 1; ++
m) {
320 auto &Xn_1 = Xn_vec[
m - 1];
321 auto &Xn_2 = Xn_vec[
n -
m];
322 diff_log_map(
i,
j,
k,
l) += cof * Xn_1(
i,
k) * Xn_2(
l,
j);
327 diff_log_map_tmp(
i,
j,
k,
l) = isddg(
i,
j,
m,
n) * diff_log_map(
m,
n,
k,
l);
335 Tensor2<T, 3, 3> &eig_vec) {
341 constexpr int lda = 3;
342 constexpr int lwork = 15;
343 std::array<double, 15> work;
344 if (
lapack_dsyev(
'V',
'U',
n, &(eigen_vec_lap(0, 0)), lda, &eig(0),
345 work.data(), lwork) > 0)
347 "The algorithm failed to compute eigenvalues.");
350 eig_vec(
i,
j) = eigen_vec_lap(
i,
j);
362template <
typename T1,
typename T2,
typename T3>
365 Tensor2<T3, 3, 3> &eig_vec) {
366 Tensor2_symmetric<double, 3> lnX;
367 Tensor2_symmetric<double, 3> diag(0., 0., 0., 0., 0., 0.);
368 for (
int ii = 0; ii != 3; ++ii)
369 diag(ii, ii) = log(
lambda(ii));
371 lnX(
i,
j) = eig_vec(
l,
i) ^ (diag(
l,
k) * eig_vec(
k,
j));
376template <
typename T>
inline auto get_ln_X(
const Tensor2_symmetric<T, 3> &X) {
380template <
typename T,
typename TP,
typename T3>
383 Tensor2_symmetric<TP, 3> &stress,
385 Tensor4<T3, 3, 3, 3, 3> &TLs3,
386 bool compute_tangent) {
389 Tensor4<double, 3, 3, 3, 3> P3;
390 Tensor2<double, 3, 3> invF;
393 CHKERR determinantTensor3by3(
F, det_f);
394 CHKERR invertTensor3by3(
F, det_f, invF);
402 MatrixDouble Vi(3, 3);
403 MatrixDouble Ksi(3, 3);
404 MatrixDouble Zeta(3, 3);
413 for (
int dd : {0, 1, 2}) {
414 for (
int cc : {0, 1, 2})
415 _N[dd](cc) = eig_vec(dd, cc);
418 for (
int dd : {0, 1, 2})
419 _n[dd](
i) =
F(
i,
j) * _N[dd](
j);
421 for (
int ii = 0; ii != 3; ++ii) {
422 e(ii) = 0.5 * log(
lambda(ii));
427 const double &lam0 =
lambda(0);
428 const double &lam1 =
lambda(1);
429 const double &lam2 =
lambda(2);
431 auto compute_ksi_eta_2eq = [&](
int i_,
int j_,
int k_) {
432 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
433 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
435 for (
int ii = 0; ii != 3; ++ii)
436 for (
int jj = 0; jj != 3; ++jj)
438 if (jj == k_ || ii == k_) {
439 Vi(ii, jj) = (e(ii) - e(jj)) / (
lambda(ii) -
lambda(jj));
449 for (
int ii = 0; ii != 3; ++ii)
450 for (
int jj = 0; jj != 3; ++jj) {
453 Vi(ii, jj) = (e(ii) - e(jj)) / (
lambda(ii) -
lambda(jj));
454 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (
lambda(ii) -
lambda(jj));
455 for (
int kk = 0; kk != 3; ++kk) {
457 if (kk != ii && kk != jj)
464 }
else if (
is_eq(lam0, lam1) &&
is_eq(lam0, lam2) &&
is_eq(lam1, lam2)) {
466 for (
int ii = 0; ii != 3; ++ii)
467 for (
int jj = 0; jj != 3; ++jj)
470 Vi(ii, jj) = 0.5 * d(ii);
471 Ksi(ii, jj) = 0.125 * f(ii);
474 }
else if (
is_eq(lam0, lam1)) {
476 compute_ksi_eta_2eq(0, 1, 2);
478 }
else if (
is_eq(lam0, lam2)) {
480 compute_ksi_eta_2eq(0, 2, 1);
482 }
else if (
is_eq(lam1, lam2)) {
484 compute_ksi_eta_2eq(1, 2, 0);
489 Tensor2<double, 3, 3> Mii;
490 Tensor2<double, 3, 3> Mij;
491 Tensor2<double, 3, 3> Mik;
492 Tensor2<double, 3, 3> Mjk;
493 Tensor2<double, 3, 3> Mjj;
495 for (
int ii = 0; ii != 3; ++ii) {
496 Mii(
i,
j) = _n[ii](
i) * _N[ii](
j) + _n[ii](
i) * _N[ii](
j);
497 P3(
i,
j,
k,
l) += 0.5 * d(ii) * _N[ii](
i) * _N[ii](
j) * Mii(
k,
l);
498 for (
int jj = 0; jj != 3; ++jj)
500 Mij(
i,
j) = _n[ii](
i) * _N[jj](
j) + _n[jj](
i) * _N[ii](
j);
501 P3(
i,
j,
k,
l) += Vi(ii, jj) * _N[ii](
i) * _N[jj](
j) * Mij(
k,
l);
505 if (!compute_tangent)
508 TLs3(
i,
j,
k,
l) = 0.;
509 for (
int ii = 0; ii != 3; ++ii)
510 for (
int jj = 0; jj != 3; ++jj)
511 Zeta(ii, jj) = (stress(
i,
j) * (_N[ii](
i) * _N[jj](
j)));
513 for (
int ii = 0; ii != 3; ++ii) {
514 Mii(
i,
j) = _n[ii](
i) * _N[ii](
j) + _n[ii](
i) * _N[ii](
j);
515 TLs3(
i,
j,
k,
l) += 0.25 * f(ii) * Zeta(ii, ii) * Mii(
i,
j) * Mii(
k,
l);
518 for (
int ii = 0; ii != 3; ++ii)
519 for (
int jj = 0; jj != 3; ++jj)
521 for (
int kk = 0; kk != 3; ++kk)
522 if (kk != ii && kk != jj) {
523 Mik(
i,
j) = _n[ii](
i) * _N[kk](
j) + _n[kk](
i) * _N[ii](
j);
524 Mjk(
i,
j) = _n[jj](
i) * _N[kk](
j) + _n[kk](
i) * _N[jj](
j);
525 TLs3(
i,
j,
k,
l) += 2. *
eta * Zeta(ii, jj) * Mik(
i,
j) * Mjk(
k,
l);
528 for (
int ii = 0; ii != 3; ++ii)
529 for (
int jj = 0; jj != 3; ++jj)
531 Mij(
i,
j) = _n[ii](
i) * _N[jj](
j) + _n[jj](
i) * _N[ii](
j);
532 Mjj(
i,
j) = _n[jj](
i) * _N[jj](
j) + _n[jj](
i) * _N[jj](
j);
533 const double k2 = 2. * Ksi(ii, jj);
534 const double zp = k2 * Zeta(ii, jj);
535 const double zp2 = k2 * Zeta(jj, jj);
536 TLs3(
i,
j,
k,
l) += zp * Mij(
i,
j) * Mjj(
k,
l) +
537 zp * Mjj(
i,
j) * Mij(
k,
l) +
538 zp2 * Mij(
i,
j) * Mij(
k,
l);
541 Tensor2<double, 3, 3> Piola;
542 Tensor2<double, 3, 3> invFPiola;
543 Tensor4<double, 3, 3, 3, 3> Ttmp3;
545 Piola(
k,
l) = nstress(
m,
n) * P3(
m,
n,
k,
l);
546 invFPiola(
i,
j) = invF(
i,
k) * Piola(
k,
j);
547 Ttmp3(
i,
j,
k,
l) = 0.5 * invFPiola(
i,
k) * kronecker_delta(
j,
l) +
548 0.5 * invFPiola(
i,
l) * kronecker_delta(
j,
k);
549 TLs3(
i,
j,
k,
l) += Ttmp3(
i,
j,
k,
l);
554template <
typename T,
typename TP,
typename T3>
558 Tensor2<T, 3, 3> &eig_vec, Ddg<T3, 3, 3> &TLs3) {
561 Tensor2<double, 3, 3> eigen_vec;
562 eigen_vec(
i,
j) = eig_vec(
i,
j);
564 auto get_uniq_nb = [&](
auto eig) {
565 return distance(&eig(0), unique(&eig(0), &eig(0) + 3,
is_eq));
568 auto sort_eigen_vals2 = [&](
auto &eig,
auto &eigen_vec) {
569 if (
is_eq(eig(0), eig(1))) {
570 Tensor2<double, 3, 3> eigen_vec_c{
571 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
572 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2),
573 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2)};
576 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
578 Tensor2<double, 3, 3> eigen_vec_c{
579 eigen_vec(1, 0), eigen_vec(1, 1), eigen_vec(1, 2),
580 eigen_vec(0, 0), eigen_vec(0, 1), eigen_vec(0, 2),
581 eigen_vec(2, 0), eigen_vec(2, 1), eigen_vec(2, 2)};
584 eigen_vec(
i,
j) = eigen_vec_c(
i,
j);
588 Tensor2_symmetric<double, 3> Ts;
589 Ts(
i,
j) = stress(
i,
j);
590 auto f = [](
double v) {
return 0.5 * log(
v); };
591 auto d_f = [](
double v) {
return 0.5 /
v; };
592 auto dd_f = [](
double v) {
return -0.5 / (
v *
v); };
594 auto nb_uniq = get_uniq_nb(eig);
596 sort_eigen_vals2(eig, eigen_vec);
605 TLs3(
i,
j,
k,
l) = t_dd(
i,
j,
k,
l);
609template <
bool TANGENT,
typename T,
typename TP,
typename T3>
612 Tensor2<T, 3, 3> &eig_vec,
613 Ddg<T3, 3, 3> &TLs3) {
615 Ddg<double, 3, 3> P3;
621 MatrixDouble Vi(3, 3);
622 MatrixDouble Ksi(3, 3);
623 MatrixDouble Zeta(3, 3);
630 for (
int dd = 0; dd != 3; ++dd) {
631 for (
int cc = 0; cc != 3; ++cc)
632 _N[dd](cc) = eig_vec(dd, cc);
635 for (
int ii = 0; ii != 3; ++ii) {
636 e(ii) = 0.5 * log(
lambda(ii));
641 const double &lam0 =
lambda(0);
642 const double &lam1 =
lambda(1);
643 const double &lam2 =
lambda(2);
645 auto compute_ksi_eta_2eq = [&](
int i_,
int j_,
int k_) {
646 Vi(i_, j_) = Vi(j_, i_) = 0.5 * d(i_);
647 Ksi(i_, j_) = Ksi(j_, i_) = 0.125 * f(i_);
649 for (
int ii = 0; ii != 3; ++ii)
650 for (
int jj = 0; jj != 3; ++jj)
652 if (jj == k_ || ii == k_) {
653 Vi(ii, jj) = (e(ii) - e(jj)) / (
lambda(ii) -
lambda(jj));
663 for (
int ii = 0; ii != 3; ++ii)
664 for (
int jj = 0; jj != 3; ++jj) {
667 Vi(ii, jj) = (e(ii) - e(jj)) / (
lambda(ii) -
lambda(jj));
668 Ksi(ii, jj) = (Vi(ii, jj) - 0.5 * d(jj)) / (
lambda(ii) -
lambda(jj));
669 for (
int kk = 0; kk != 3; ++kk) {
671 if (kk != ii && kk != jj)
678 }
else if (
is_eq(lam0, lam1) &&
is_eq(lam0, lam2) &&
is_eq(lam1, lam2)) {
680 for (
int ii = 0; ii != 3; ++ii)
681 for (
int jj = 0; jj != 3; ++jj)
684 Vi(ii, jj) = 0.5 * d(ii);
685 Ksi(ii, jj) = 0.125 * f(ii);
688 }
else if (
is_eq(lam0, lam1)) {
690 compute_ksi_eta_2eq(0, 1, 2);
692 }
else if (
is_eq(lam0, lam2)) {
694 compute_ksi_eta_2eq(0, 2, 1);
696 }
else if (
is_eq(lam1, lam2)) {
698 compute_ksi_eta_2eq(1, 2, 0);
702 Tensor2_symmetric<double, 3> Mii;
703 Tensor2_symmetric<double, 3> Mij;
704 Tensor2_symmetric<double, 3> Mik;
705 Tensor2_symmetric<double, 3> Mjk;
706 Tensor2_symmetric<double, 3> Mjj;
708 for (
int ii = 0; ii != 3; ++ii) {
709 Mii(
i,
j) = (_N[ii](
i) ^ _N[ii](
j)) + (_N[ii](
i) ^ _N[ii](
j));
710 P3(
i,
j,
k,
l) += ((0.5 * d(ii) * _N[ii](
i)) ^ _N[ii](
j)) * Mii(
k,
l);
711 for (
int jj = 0; jj != 3; ++jj)
713 Mij(
i,
j) = (_N[ii](
i) ^ _N[jj](
j)) + (_N[jj](
i) ^ _N[ii](
j));
714 P3(
i,
j,
k,
l) += ((Vi(ii, jj) * _N[ii](
i)) ^ _N[jj](
j)) * Mij(
k,
l);
721 TLs3(
i,
j,
k,
l) = 0.;
722 for (
int ii = 0; ii != 3; ++ii)
723 for (
int jj = 0; jj != 3; ++jj)
724 Zeta(ii, jj) = stress(
i,
j) * (_N[ii](
i) ^ _N[jj](
j));
726 for (
int ii = 0; ii != 3; ++ii) {
727 Mii(
i,
j) = (_N[ii](
i) ^ _N[ii](
j)) + (_N[ii](
i) ^ _N[ii](
j));
728 TLs3(
i,
j,
k,
l) += (0.25 * f(ii) * Zeta(ii, ii) * Mii(
i,
j)) * Mii(
k,
l);
731 for (
int ii = 0; ii != 3; ++ii)
732 for (
int jj = 0; jj != 3; ++jj)
734 Mij(
i,
j) = (_N[ii](
i) ^ _N[jj](
j)) + (_N[jj](
i) ^ _N[ii](
j));
735 Mjj(
i,
j) = (_N[jj](
i) ^ _N[jj](
j)) + (_N[jj](
i) ^ _N[jj](
j));
736 const double k2 = 2. * Ksi(ii, jj);
737 const double zp = k2 * Zeta(ii, jj);
738 const double zp2 = k2 * Zeta(jj, jj);
739 TLs3(
i,
j,
k,
l) += (zp * Mij(
i,
j)) * Mjj(
k,
l) +
740 (zp * Mjj(
i,
j)) * Mij(
k,
l) +
741 (zp2 * Mij(
i,
j)) * Mij(
k,
l);
742 for (
int kk = 0; kk != 3; ++kk) {
743 if (kk != ii && kk != jj) {
744 Mik(
i,
j) = (_N[ii](
i) ^ _N[kk](
j)) + (_N[kk](
i) ^ _N[ii](
j));
745 Mjk(
i,
j) = (_N[jj](
i) ^ _N[kk](
j)) + (_N[kk](
i) ^ _N[jj](
j));
747 (2. *
eta * Zeta(ii, jj) * Mik(
i,
j)) * Mjk(
k,
l);
756inline Tensor2<double, 3, 3>
758 Tensor2<T, 3, 3> &
F, Tensor2<double, 3, 3> &t_stress) {
760 Tensor2<double, 3, 3> out;
762 Tensor2<double, 3, 3> strainPl;
763 Tensor2<double, 3, 3> strainMin;
764 strainPl(
i,
j) =
F(
i,
j);
765 strainMin(
i,
j) =
F(
i,
j);
766 const double h = 1e-8;
767 strainPl(kk, ll) +=
h;
768 strainMin(kk, ll) -=
h;
769 Tensor4<PackPtr<double *, 1>, 3, 3, 3, 3> dummy;
771 auto calculate_stress = [&](
auto &
F) {
772 Tensor2<double, 3, 3> Piola;
773 Tensor2_symmetric<double, 3> stress_sym;
774 Tensor2_symmetric<double, 3> C;
775 Tensor2_symmetric<double, 3> strain;
780 Tensor2<double, 3, 3> eig_vec;
785 strain(
i,
j) = 0.5 * lnC(
i,
j);
786 auto t_D = getFTensor4DdgFromMat<3, 3, 0>(*common_data.
mtD);
787 stress_sym(
i,
j) = t_D(
i,
j,
k,
l) * strain(
k,
l);
801 Piola(
k,
l) = stress_tmp(
i,
j) * D2(
i,
j,
k,
l);
806 auto StressPlusH = calculate_stress(strainPl);
807 auto StressMinusH = calculate_stress(strainMin);
810 out(
i,
j) = (StressPlusH(
i,
j) - StressMinusH(
i,
j)) / (2 *
h);
816inline Tensor4<double, 3, 3, 3, 3>
818 Tensor2<double, 3, 3> &t_stress) {
819 Tensor4<double, 3, 3, 3, 3> my_D;
820 for (
int k = 0;
k != 3; ++
k) {
821 for (
int l = 0;
l != 3; ++
l) {
823 for (
int i = 0;
i != 3; ++
i) {
824 for (
int j = 0;
j != 3; ++
j) {
825 my_D(
i,
j,
k,
l) = d_stress(
i,
j);
834typedef boost::function<Tensor1<double, 3>(
const double,
const double,
840 boost::shared_ptr<CommonData> common_data_ptr);
849 boost::shared_ptr<CommonData> common_data_ptr);
858 boost::shared_ptr<CommonData> common_data_ptr, boost::shared_ptr<MatrixDouble> mat_D_Ptr);
866template <
bool TANGENT>
869 boost::shared_ptr<CommonData> common_data_ptr,
870 boost::shared_ptr<MatrixDouble> mat_D_ptr);
881 boost::shared_ptr<CommonData> common_data_ptr);
888template <
bool LOGSTRAIN,
bool REACTIONS>
891 boost::shared_ptr<CommonData> common_data_ptr);
901 moab::Interface &post_proc_mesh,
902 std::vector<EntityHandle> &map_gauss_pts,
903 boost::shared_ptr<CommonData> common_data_ptr);
914 boost::shared_ptr<CommonData> common_data_ptr,
915 std::map<EntityHandle, int> &map_block_tets);
926 boost::shared_ptr<CommonData> common_data_ptr,
927 std::map<EntityHandle, int> &map_block_tets);
937 moab::Interface &post_proc_mesh,
938 std::vector<EntityHandle> &map_gauss_pts,
939 boost::shared_ptr<CommonData> common_data_ptr);
952 boost::shared_ptr<CommonData> common_data_ptr,
953 VectorDouble &disp,
const Range &essential_ents)
958 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &row_data) {
960 const int nb_row_dofs = row_data.getIndices().size();
961 const int nb_row_base_functions = row_data.getN().size2();
970 auto t_disp = getFTensor1FromMat<3>(*
commonDataPtr->mDispPtr);
972 const int nb_integration_pts = getGaussPts().size2();
973 auto t_row_val = row_data.getFTensor0N();
974 auto t_w = getFTensor0IntegrationWeight();
975 const double vol = getMeasure();
976 auto t_coords = getFTensor1CoordsAtGaussPts();
979 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
980 const double a = vol * t_w;
987 auto t_nf = getFTensor1FromArray<3, 3>(locF);
989 for (;rr != nb_row_dofs / 3; ++rr) {
990 t_nf(
i) +=
a * t_row_val * (t_disp(
i) - bc_disp(
i));
994 for (; rr < nb_row_base_functions; ++rr)
1013 boost::shared_ptr<VectorDouble> c_coords,
1014 boost::shared_ptr<Range> ents_ptr)
1018 MoFEMErrorCode
iNtegrate(EntitiesFieldData::EntData &row_data) {
1020 const int nb_row_dofs = row_data.getIndices().size();
1021 const int nb_row_base_functions = row_data.getN().size2();
1030 auto t_coords = getFTensor1CoordsAtGaussPts();
1032 const int nb_integration_pts = getGaussPts().size2();
1033 auto t_row_val = row_data.getFTensor0N();
1034 auto t_w = getFTensor0IntegrationWeight();
1035 const double vol = getMeasure();
1037 auto t_omega = getFTensor1FromArray<3, 3>(
sourceVec);
1038 auto t_cen_coords = getFTensor1FromArray<3, 3>(*
centerCoords);
1043 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1044 const double a = vol * t_w;
1045 c_dist(
i) = t_cen_coords(
i) - t_coords(
i);
1046 rot_disp(
i) = c_dist(
i) - rot(
j,
i) * c_dist(
j);
1048 auto t_nf = getFTensor1FromArray<3, 3>(locF);
1050 for (; rr != nb_row_dofs / 3; ++rr) {
1051 t_nf(
i) -=
a * t_row_val * rot_disp(
i);
1055 for (; rr < nb_row_base_functions; ++rr)
1075 boost::shared_ptr<CommonData> common_data_ptr)
1084 if (
auto stored_data_ptr =
1085 e_ptr->getSharedStoragePtr<EssentialBcStorage>()) {
1088 data.
getIndices() = stored_data_ptr->entityIndices;
1104 boost::shared_ptr<CommonData> common_data_ptr)
1115 boost::shared_ptr<CommonData> common_data_ptr,
1116 const double eps = 1e-12)
1139 boost::shared_ptr<CommonData> common_data_ptr)
1150 boost::shared_ptr<CommonData> common_data_ptr,
1151 const double eps = 1e-12)
1172template <
typename TYPE>
1175 const std::string col_field_name, ScalarFun beta,
1176 boost::shared_ptr<CommonData> common_data_ptr)
1177 : TYPE(row_field_name, col_field_name, beta),
1199 const std::string col_field_name, ScalarFun beta,
1200 boost::shared_ptr<CommonData> common_data_ptr)
1236 boost::ptr_vector<DataFromMove> &bc_data,
1237 boost::shared_ptr<CommonData> common_data_ptr)
1243 CHKERR VecSetOption(ts_F, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1246 case CTX_TSSETIJACOBIAN: {
1247 snes_ctx = CTX_SNESSETJACOBIAN;
1259 case CTX_TSSETIFUNCTION: {
1260 snes_ctx = CTX_SNESSETFUNCTION;
1273 for (
auto &bdata :
bcData) {
1274 bdata.scaledValues = bdata.values;
1276 bdata.scaledValues);
1280 if (time <= (*cache).rollers_stop_time) {
1283 for (
auto &roller : (*cache).rollerDataVec) {
1284 roller.BodyDispScaled = roller.rollerDisp;
1289 for (
auto &roller : (*cache).rollerDataVec) {
1290 if (roller.methodOpForRollerPosition) {
1291 roller.BodyDispScaled.clear();
1293 this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
1295 if (roller.methodOpForRollerDirection) {
1296 roller.BodyDirectionScaled.clear();
1298 this, roller.methodOpForRollerDirection,
1299 roller.BodyDirectionScaled);
1321 auto assemble = [](Mat
a) {
1324 CHKERR MatAssemblyBegin(
a, MAT_FINAL_ASSEMBLY);
1325 CHKERR MatAssemblyEnd(
a, MAT_FINAL_ASSEMBLY);
1330 if (
ts_ctx == CTX_TSSETIJACOBIAN) {
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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 ...
#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 ...
constexpr double omega
Save field DOFS on vertices/tags.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
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::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)
constexpr double t
plate stiffness
constexpr auto field_name
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorInt & getIndices() const
Get global indices of dofs on entity.
default operator for TRI element
@ OPROW
operator doWork function is executed on FE rows
Not used at this stage. Could be used to do some calculations, before assembly of local elements.
MoFEMErrorCode preProcess()
boost::shared_ptr< MethodForForceScaling > methodsOpForMove
boost::ptr_vector< DataFromMove > & bcData
MoFEMErrorCode postProcess()
FePrePostProcess(MoFEM::Interface &m_field, boost::ptr_vector< DataFromMove > &bc_data, boost::shared_ptr< CommonData > common_data_ptr)
MoFEM::Interface & mField
boost::shared_ptr< MethodForForceScaling > methodsOpForRollersDisp
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
std::vector< EntityHandle > & mapGaussPts
moab::Interface & postProcMesh
boost::shared_ptr< CommonData > commonDataPtr