386 {
387
388
389 Tensor4<double, 3, 3, 3, 3> P3;
390 Tensor2<double, 3, 3> invF;
391
392 double det_f;
393 CHKERR determinantTensor3by3(
F, det_f);
394 CHKERR invertTensor3by3(
F, det_f, invF);
395
396
397
398
399 VectorDouble e(3);
400 VectorDouble f(3);
401 VectorDouble d(3);
402 MatrixDouble Vi(3, 3);
403 MatrixDouble Ksi(3, 3);
404 MatrixDouble Zeta(3, 3);
405 Vi.clear();
406 Ksi.clear();
407 Zeta.clear();
409
412
413 for (int dd : {0, 1, 2}) {
414 for (int cc : {0, 1, 2})
415 _N[dd](cc) = eig_vec(dd, cc);
416 }
417
418 for (int dd : {0, 1, 2})
419 _n[dd](
i) =
F(
i,
j) * _N[
dd](
j);
420
421 for (int ii = 0; ii != 3; ++ii) {
422 e(ii) = 0.5 * log(
lambda(ii));
425 }
426
427 const double &lam0 =
lambda(0);
428 const double &lam1 =
lambda(1);
429 const double &lam2 =
lambda(2);
430
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_);
434
435 for (int ii = 0; ii != 3; ++ii)
436 for (int jj = 0; jj != 3; ++jj)
437 if (ii != jj) {
438 if (jj == k_ || ii == k_) {
439 Vi(ii, jj) = (e(ii) - e(jj)) / (
lambda(ii) -
lambda(jj));
440 Ksi(ii, jj) =
442 }
443 }
445 };
446
448
449 for (int ii = 0; ii != 3; ++ii)
450 for (int jj = 0; jj != 3; ++jj) {
451 if (ii != jj) {
452
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) {
456
457 if (kk != ii && kk != jj)
460 }
461 }
462 }
463
464 }
else if (
is_eq(lam0, lam1) &&
is_eq(lam0, lam2) &&
is_eq(lam1, lam2)) {
465
466 for (int ii = 0; ii != 3; ++ii)
467 for (int jj = 0; jj != 3; ++jj)
468 if (ii != jj) {
469
470 Vi(ii, jj) = 0.5 *
d(ii);
471 Ksi(ii, jj) = 0.125 *
f(ii);
472 }
474 }
else if (
is_eq(lam0, lam1)) {
475
476 compute_ksi_eta_2eq(0, 1, 2);
477
478 }
else if (
is_eq(lam0, lam2)) {
479
480 compute_ksi_eta_2eq(0, 2, 1);
481
482 }
else if (
is_eq(lam1, lam2)) {
483
484 compute_ksi_eta_2eq(1, 2, 0);
485 }
486
488
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;
494
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)
499 if (ii != 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);
502 }
503 }
504
505 if (!compute_tangent)
506 return P3;
507
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)));
512
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);
516 }
517
518 for (int ii = 0; ii != 3; ++ii)
519 for (int jj = 0; jj != 3; ++jj)
520 if (jj != ii)
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);
526 }
527
528 for (int ii = 0; ii != 3; ++ii)
529 for (int jj = 0; jj != 3; ++jj)
530 if (jj != ii) {
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);
539 }
540
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);
549 TLs3(
i,
j,
k,
l) += Ttmp3(
i,
j,
k,
l);
550
551 return P3;
552};
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)