Operator for linear form, usually to calculate values on right hand side.
147 {
149
152
153
163
164
166
167 auto t_L = FTensor::SymmLTensor<double, 3>();
168
170 *
dataAtPts->getStretchTensorAtPts(), nb_integration_pts);
172 *
dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
174 *
dataAtPts->getStretchH1AtPts(), nb_integration_pts);
176 *
dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
178 *
dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
180 *
dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
182 *
dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
184 *
dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
185
187 *
dataAtPts->getDeformationGradient(), nb_integration_pts);
189 *
dataAtPts->getPlasticH(), nb_integration_pts);
191 *
dataAtPts->getPlasticF(), nb_integration_pts);
193 *
dataAtPts->getInvPlasticF(), nb_integration_pts);
195 dataAtPts->hdOmegaAtPts, nb_integration_pts);
197 dataAtPts->hdLogStretchAtPts, nb_integration_pts);
198
200 dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
202 dataAtPts->leviKirchhoff0AtPts, nb_integration_pts);
204 dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
206 dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
208 dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
209
211 dataAtPts->rotMatAtPts, nb_integration_pts);
213 *
dataAtPts->getEigenVals(), nb_integration_pts);
215 *
dataAtPts->getEigenVecs(), nb_integration_pts);
216 dataAtPts->nbUniq.resize(nb_integration_pts,
false);
218 dataAtPts->eigenValsC, nb_integration_pts);
220 dataAtPts->eigenVecsC, nb_integration_pts);
221 dataAtPts->nbUniqC.resize(nb_integration_pts,
false);
222
224 dataAtPts->logStretch2H1AtPts, nb_integration_pts);
226 dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
227
229 dataAtPts->internalStressAtPts, nb_integration_pts);
231
232
233
234
235 auto t_log_plasticH =
dataAtPts->getFTensorPlasticH(nb_integration_pts);
236 auto t_plasticF_reconstruct =
237 dataAtPts->getFTensorPlasticF(nb_integration_pts);
238 auto t_invPlasticF_reconstruct =
239 dataAtPts->getFTensorInvPlasticF(nb_integration_pts);
242 };
245 };
246 for (int gg = 0; gg != nb_integration_pts; ++gg) {
249 t_eigen_vecs(
i,
j) = t_log_plasticH(
i,
j);
252 "Failed to diagonalise logarithmic plastic deformation");
253 const auto t_exp_plasticH =
255 const auto t_inv_exp_plasticH =
257 t_plasticF_reconstruct(
i,
j) = t_exp_plasticH(
i,
j);
258 t_invPlasticF_reconstruct(
i,
j) = t_inv_exp_plasticH(
i,
j);
259
260#ifndef NDEBUG
261 const double det_plasticF =
263 if (!std::isfinite(det_plasticF) ||
264 det_plasticF <= std::numeric_limits<double>::epsilon())
266 "Plastic deformation gradient must have a positive determinant; "
267 "got %g",
268 det_plasticF);
269#endif
270
271 ++t_log_plasticH;
272 ++t_plasticF_reconstruct;
273 ++t_invPlasticF_reconstruct;
274 }
275
278 auto get_intermediate_p =
280 DL>::size(intermediate_p_at_pts, nb_integration_pts);
281 auto get_intermediate_p0 =
283 DL>::size(intermediate_p0_at_pts, nb_integration_pts);
284 auto t_reference_P =
dataAtPts->getFTensorApproxP(nb_integration_pts);
285 auto t_reference_P0 =
dataAtPts->getFTensorApproxP0(nb_integration_pts);
286 auto t_plasticF =
dataAtPts->getFTensorPlasticF(nb_integration_pts);
287 auto t_intermediate_P = get_intermediate_p();
288 auto t_intermediate_P0 = get_intermediate_p0();
289 for (int gg = 0; gg != nb_integration_pts; ++gg) {
291 t_intermediate_P(
i,
j) =
292 t_reference_P(
i,
k) * t_plasticF(
j,
k) / det_plasticF;
293 t_intermediate_P0(
i,
j) =
294 t_reference_P0(
i,
k) * t_plasticF(
j,
k) / det_plasticF;
295 ++t_reference_P;
296 ++t_reference_P0;
297 ++t_plasticF;
298 ++t_intermediate_P;
299 ++t_intermediate_P0;
300 }
301
302
305 auto t_h_dlog_u =
307 auto t_levi_kirchhoff =
309 auto t_levi_kirchhoff0 =
311 auto t_levi_kirchhoff_domega =
313 auto t_levi_kirchhoff_dstreach =
315 auto t_levi_kirchhoff_dP =
317 auto t_approx_P_adjoint_dstretch =
319 auto t_approx_P_adjoint_log_du =
321 auto t_approx_P_adjoint_log_du_dP =
323 auto t_approx_P_adjoint_log_du_domega =
331 auto t_nb_uniq =
333 auto t_eigen_vals_C =
dataAtPts->getFTensorEigenValsC(nb_integration_pts);
334 auto t_eigen_vecs_C =
dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
336 auto t_nb_uniq_C =
338
341 auto t_log_stretch_total =
344
345
348 auto t_approx_P = get_intermediate_p();
349 auto t_approx_P0 = get_intermediate_p0();
351
352
355
356 auto next = [&]() {
357
358 ++t_h;
359 ++t_h_domega;
360 ++t_h_dlog_u;
361 ++t_levi_kirchhoff;
362 ++t_levi_kirchhoff0;
363 ++t_levi_kirchhoff_domega;
364 ++t_levi_kirchhoff_dstreach;
365 ++t_levi_kirchhoff_dP;
366 ++t_approx_P_adjoint_dstretch;
367 ++t_approx_P_adjoint_log_du;
368 ++t_approx_P_adjoint_log_du_dP;
369 ++t_approx_P_adjoint_log_du_domega;
370 ++t_R;
371 ++t_u;
372 ++t_diff_u;
373 ++t_eigen_vals;
374 ++t_eigen_vecs;
375 ++t_nb_uniq;
376 ++t_eigen_vals_C;
377 ++t_eigen_vecs_C;
378 ++t_nb_uniq_C;
379 ++t_u_h1;
380 ++t_diff_u_h1;
381 ++t_log_u2_h1;
382 ++t_log_stretch_total;
383
384 ++t_omega;
385 ++t_omega0;
386 ++t_grad_h1;
387 ++t_approx_P;
388 ++t_approx_P0;
389 ++t_log_u;
390 ++t_log_u0;
391 };
392
395 constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
396
397 auto calculate_stretch_from_log = [&](auto &t_log_u_src, auto &t_u_dst,
398 auto &t_eigen_vals_dst,
399 auto &t_eigen_vecs_dst,
400 int &nb_uniq_dst) {
404 eigen_vec(
i,
j) = t_log_u_src(
i,
j);
406 MOFEM_LOG(
"SELF", Sev::error) <<
"Failed to compute eigen values";
407 }
408
409
410 nb_uniq_dst = getUniqNb<3>(eig);
411 if (nb_uniq_dst < 3) {
412 CHKERR sortEigenVals<3>(eig, eigen_vec);
413 }
414 t_eigen_vals_dst(
i) = eig(
i);
415 t_eigen_vecs_dst(
i,
j) = eigen_vec(
i,
j);
419 };
420
421 auto calculate_log_stretch = [&]() {
423 int nb_uniq_val = 0;
424 CHKERR calculate_stretch_from_log(t_log_u, t_u, t_eigen_vals, t_eigen_vecs,
425 nb_uniq_val);
426 t_nb_uniq = nb_uniq_val;
427 auto get_t_diff_u = [&]() {
430 t_nb_uniq);
431 };
432 t_diff_u(
i,
j,
k,
l) = get_t_diff_u()(
i,
j,
k,
l);
434 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
436 };
437
438 auto calculate_total_stretch = [&](auto &t_h1) {
441
442 t_log_u2_h1(
i,
j) = 0;
443 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
444
445 } else {
446
450
452 t_C_h1(
i,
j) = t_h1(
k,
i) * t_h1(
k,
j);
453 t_eigen_vec(
i,
j) = t_C_h1(
i,
j);
456 "Failed to compute eigenvalues of F_H1^T F_H1");
457 }
458
459 t_nb_uniq_C = getUniqNb<3>(t_eig_C);
460 if (t_nb_uniq_C < 3) {
461 CHKERR sortEigenVals<3>(t_eig_C, t_eigen_vec);
462 }
463 for (int aa = 0; aa != 3; ++aa) {
464 if (!std::isfinite(t_eig_C(aa)) || t_eig_C(aa) <= 0.) {
466 "F_H1^T F_H1 must be positive definite; eigenvalue %d is "
467 "%g",
468 aa, t_eig_C(aa));
469 }
470 const double principal_stretch = std::sqrt(t_eig_C(aa));
471 const double coordinate_stretch =
473 if (!std::isfinite(coordinate_stretch)) {
474 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
475 "Non-finite H1 coordinate stretch for principal stretch %g",
476 principal_stretch);
477 }
478 t_coordinate_stretch(aa) = coordinate_stretch;
479 }
480 t_eigen_vals_C(
i) = t_eig_C(
i);
481 t_eigen_vecs_C(
i,
j) = t_eigen_vec(
i,
j);
482
485 [](
const double v) {
return v; })(
i,
j);
486
487
488 t_log_stretch_total(
i,
j) = t_log_u2_h1(
i,
j) + t_log_u(
i,
j);
489 }
491 };
492
493 auto no_h1_loop = [&]() {
495
498 break;
500 break;
501 default:
503 "no_h1_loop is only implemented for LARGE_ROT");
504 };
505
506 for (int gg = 0; gg != nb_integration_pts; ++gg) {
507
509
512
513
514 CHKERR calculate_log_stretch();
516 if ((
dataAtPts->physicsPtr->getFeatures() &
517 PhysicalEquations::noStretchMask)
518 .any()) {
519 t_u0(
i,
j) = t_u(
i,
j);
520 } else {
523 int nb_uniq_0 = 0;
524 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
525 t_eigen_vecs_0, nb_uniq_0);
526 }
527
528 CHKERR calculate_total_stretch(t_h1);
529
530 t_u_h1(
i,
j) = t_u(
i,
j);
531 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i,
j,
k,
l);
533 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
534
537
538 auto large_rot = [&]() {
542 t_diff_diff_R(
i,
j,
k,
l) =
544
548
549 t_h(
i,
k) = t_R(
i,
l) * t_u(
l,
k);
550
552 t_rotated_P(
l,
k) = t_R(
i,
l) * t_approx_P(
i,
k);
553 t_approx_P_adjoint_dstretch(
l,
k) =
554 t_diff_sym(
l,
k,
i,
j) * t_rotated_P(
i,
j);
555 t_approx_P_adjoint_log_du(
L) =
556 t_approx_P_adjoint_dstretch(
l,
k) * t_Ldiff_u(
l,
k,
L);
557
558 t_levi_kirchhoff(
m) =
559 t_diff_R(
i,
l,
m) * (t_u(
l,
k) * t_approx_P(
i,
k));
560 t_levi_kirchhoff0(
m) =
561 t_diff_R0(
i,
l,
m) * (t_u0(
l,
k) * t_approx_P0(
i,
k));
562
564 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
565 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
566
567 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
568 t_R(
i,
l) * t_Ldiff_u(
l,
k,
L);
569
571 t_A(
k,
l,
m) = t_diff_R(
i,
l,
m) * t_approx_P(
i,
k);
572 t_approx_P_adjoint_log_du_domega(
m,
L) =
573 t_A(
k,
l,
m) * t_Ldiff_u(
k,
l,
L);
574
575 t_levi_kirchhoff_dstreach(
m,
L) =
576 t_diff_R(
i,
l,
m) * (t_Ldiff_u(
l,
k,
L) * t_approx_P(
i,
k));
577 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u(
l,
k);
578 t_levi_kirchhoff_domega(
m,
n) =
579 t_diff_diff_R(
i,
l,
m,
n) * (t_u(
l,
k) * t_approx_P(
i,
k));
580
581 if (
dataAtPts->physicsPtr->getFeatures().test(
582 PhysicalEquations::NO_STRETCH_NONLINEAR)) {
586 auto [t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega] =
588 t_approx_P);
589
590 t_h_domega(
i,
k,
m) += t_d_h_d_omega(
i,
k,
m);
592 t_d_u_contract_p;
593 t_d_u_contract_p(
i,
l,
n) =
594 t_d_u_d_omega(
l,
k,
n) * t_approx_P(
i,
k);
595 t_levi_kirchhoff_domega(
m,
n) +=
596 t_diff_R(
i,
l,
m) * t_d_u_contract_p(
i,
l,
n);
597
601 t_d_b_d_p;
602 t_d_b_d_p(
i,
j,
k,
l) =
603 t_diff_sym(
i,
j,
m,
l) * t_R(
k,
m);
606 t_d_u_d_p;
607 t_d_u_d_p(
i,
j,
k,
l) =
608 t_d_u_d_b(
i,
j,
m,
n) * t_d_b_d_p(
m,
n,
k,
l);
609 t_levi_kirchhoff_dP(
m,
k,
l) +=
610 t_d_u_d_p(
i,
j,
k,
l) * t_d_b_d_omega(
i,
j,
m);
611 }
612 }
613 }
614 };
615
616 auto moderate_rot = [&](auto &t_omega0) {
618 "moderate_rot is not implemented yet");
619 };
620
621 auto small_rot = [&]() {
622 t_u_h1(
i,
j) = t_u(
i,
j);
623 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i,
j,
k,
l);
625 t_Ldiff_u(
i,
j,
L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n,
L);
626
629
631 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
632
633
635 t_rotated_P(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
636 t_approx_P_adjoint_dstretch(
i,
j) =
637 t_diff_sym(
i,
j,
k,
l) * t_rotated_P(
k,
l);
638 t_approx_P_adjoint_log_du(
L) =
639 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
640 t_approx_P_adjoint_log_du_dP(
i,
j,
L) =
641 t_R(
i,
k) * t_Ldiff_u(
k,
j,
L);
642 t_approx_P_adjoint_log_du_domega(
m,
L) =
645
646
649 t_levi_kirchhoff_dstreach(
m,
L) = 0;
651 t_levi_kirchhoff_domega(
m,
n) = 0;
652 };
653
654
657 large_rot();
658 break;
660 moderate_rot(t_omega0);
661 break;
663 small_rot();
664 break;
665 default:
667 "rotationSelector not handled");
668 }
669
670 next();
671 }
672
674 };
675
676 auto large_loop = [&]() {
678
681 break;
683 break;
684 default:
686 "rotSelector should be large or small");
687 };
688
689 for (int gg = 0; gg != nb_integration_pts; ++gg) {
690
692
697 break;
698 default:
700 "Selected grad approximator not handled");
701 };
702
703
704 CHKERR calculate_log_stretch();
706 if ((
dataAtPts->physicsPtr->getFeatures() &
707 PhysicalEquations::noStretchMask)
708 .any()) {
709 t_u0(
i,
j) = t_u(
i,
j);
710 } else {
713 int nb_uniq_0 = 0;
714 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
715 t_eigen_vecs_0, nb_uniq_0);
716 }
717
718 CHKERR calculate_total_stretch(t_h1);
719
720 t_u_h1(
l,
k) = t_u(
l, o) * t_h1(o,
k);
722 t_u_h10(
l,
k) = t_u0(
l, o) * t_h1(o,
k);
723 t_diff_u_h1(
i,
j,
k,
l) = t_diff_u(
i, o,
k,
l) * t_h1(o,
j);
725 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
726
730
731
737 t_diff_diff_R(
i,
j,
l,
m) = 0;
738 break;
745 t_diff_diff_R(
i,
j,
k,
l) =
747 break;
748
749 default:
751 "rotationSelector not handled");
752 }
753
754
755 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
756
757
760 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
761 t_approx_P_adjoint_dstretch(
l, o) =
762 t_diff_sym(
l, o,
i,
j) * t_rotated_P(
i,
j);
763 t_approx_P_adjoint_log_du(
L) =
764 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
765
766
767 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
768 t_levi_kirchhoff0(
m) =
769 t_diff_R0(
i,
l,
m) * t_u_h10(
l,
k) * t_approx_P0(
i,
k);
770
772
773 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
774 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
775
776 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
777 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
778
780 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
781 t_approx_P_adjoint_log_du_domega(
m,
L) =
782 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
783
784 t_levi_kirchhoff_dstreach(
m,
L) =
785 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
786
787 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
788 t_levi_kirchhoff_domega(
m,
n) =
789 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
790 }
791
792 next();
793 }
794
796 };
797
798 auto moderate_loop = [&]() {
800
803 break;
805 break;
806 default:
808 "rotSelector should be large or small");
809 };
810
811 for (int gg = 0; gg != nb_integration_pts; ++gg) {
812
814
819 break;
820 default:
822 "Selected grad approximator not handled");
823 };
824
825
826 CHKERR calculate_log_stretch();
827
828 CHKERR calculate_total_stretch(t_h1);
829
830 auto t_diff = FTensor::DiffTensor<double>();
831
832 t_u_h1(
l,
k) = (
t_kd(
l, o) + t_log_u(
l, o)) * t_h1(o,
k);
834 t_u_h10(
l,
k) = (
t_kd(
l, o) + t_log_u0(
l, o)) * t_h1(o,
k);
835 t_diff_u_h1(
i,
j,
k,
l) = t_diff(
i, o,
k,
l) * t_h1(o,
j);
837 t_Ldiff_u_h1(
l,
k,
L) = t_diff_u_h1(
l,
k,
i,
j) * t_L(
i,
j,
L);
838
842
843
849 t_diff_diff_R(
i,
j,
l,
m) = 0;
850 break;
857 t_diff_diff_R(
i,
j,
k,
l) =
859 break;
860
861 default:
863 "rotationSelector not handled");
864 }
865
866
867 t_h(
i,
k) = t_R(
i,
l) * t_u_h1(
l,
k);
868
869
872 (t_R(
i,
l) * t_approx_P(
i,
k)) * t_h1(o,
k);
873 t_approx_P_adjoint_dstretch(
l, o) =
874 t_diff_sym(
l, o,
i,
j) * t_rotated_P(
i,
j);
875 t_approx_P_adjoint_log_du(
L) =
876 t_R(
i,
l) * t_approx_P(
i,
k) * t_Ldiff_u_h1(
l,
k,
L);
877
878
879 t_levi_kirchhoff(
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k) * t_approx_P(
i,
k);
880 t_levi_kirchhoff0(
m) =
881 t_diff_R0(
i,
l,
m) * t_u_h10(
l,
k) * t_approx_P0(
i,
k);
882
884
885 t_h_domega(
i,
k,
m) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
886 t_h_dlog_u(
i,
k,
L) = t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
887
888 t_approx_P_adjoint_log_du_dP(
i,
k,
L) =
889 t_R(
i,
l) * t_Ldiff_u_h1(
l,
k,
L);
890
892 t_A(
m,
L,
i,
k) = t_diff_R(
i,
l,
m) * t_Ldiff_u_h1(
l,
k,
L);
893 t_approx_P_adjoint_log_du_domega(
m,
L) =
894 t_A(
m,
L,
i,
k) * t_approx_P(
i,
k);
895
896 t_levi_kirchhoff_dstreach(
m,
L) =
897 t_diff_R(
i,
l,
m) * (t_Ldiff_u_h1(
l,
k,
L) * t_approx_P(
i,
k));
898
899 t_levi_kirchhoff_dP(
m,
i,
k) = t_diff_R(
i,
l,
m) * t_u_h1(
l,
k);
900 t_levi_kirchhoff_domega(
m,
n) =
901 t_diff_diff_R(
i,
l,
m,
n) * (t_u_h1(
l,
k) * t_approx_P(
i,
k));
902 }
903
904 next();
905 }
906
908 };
909
910 auto small_loop = [&]() {
914 break;
915 default:
917 "rotSelector should be small");
918 };
919
920 for (int gg = 0; gg != nb_integration_pts; ++gg) {
921
926 break;
927 default:
929 "gradApproximator not handled");
930 };
931
935 "stretchSelector should be linear for small loop");
936 } else {
937 t_u(
i,
j) = t_symm_kd(
i,
j) + t_log_u(
i,
j);
938 t_u_h1(
i,
j) = t_u(
i,
j);
939 t_diff_u_h1(
i,
j,
k,
l) =
941 t_diff_u_h1(
i,
j,
k,
l) /= 2.;
942 t_Ldiff_u(
i,
j,
L) = t_L(
i,
j,
L);
943 }
944 t_log_u2_h1(
i,
j) = 0;
945 t_log_stretch_total(
i,
j) = t_log_u(
i,
j);
946
949
951 t_h_dlog_u(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
952
953
954 t_approx_P_adjoint_dstretch(
i,
j) =
955 t_diff_sym(
i,
j,
k,
l) * t_approx_P(
k,
l);
956 t_approx_P_adjoint_log_du(
L) =
957 t_approx_P_adjoint_dstretch(
i,
j) * t_Ldiff_u(
i,
j,
L);
958 t_approx_P_adjoint_log_du_dP(
i,
j,
L) = t_Ldiff_u(
i,
j,
L);
959 t_approx_P_adjoint_log_du_domega(
m,
L) = 0;
960
961
964 t_levi_kirchhoff_dstreach(
m,
L) = 0;
966 t_levi_kirchhoff_domega(
m,
n) = 0;
967
968 next();
969 }
970
972 };
973
977 break;
980 break;
983 break;
986 break;
987 default:
989 "gradApproximator not handled");
990 break;
991 };
992
994}
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#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 MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
boost::function< T(const T)> Fun
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
auto getDiffSpatialGradientDR(TInvD &t_d_u_d_b, TRotation &t_R, TDiffRotation &t_diff_R, TStress &t_P)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
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.
decltype(GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
constexpr IntegrationType I
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static enum RotSelector rotSelector
static enum RotSelector gradApproximator
static constexpr enum SymmetrySelector symmetrySelector
static boost::function< double(const double)> f
static boost::function< double(const double)> d_f
static boost::function< double(const double)> inv_f
static auto diffDiffExp(A &&t_w_vee, B &&theta)
static auto diffExp(A &&t_w_vee, B &&theta)
static auto exp(A &&t_w_vee, B &&theta)
const TSMethod::TSContext getTSCtx() const
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
Public Member Functions inherited from MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator