324 assert(data.
x.cols() == 1);
326 const int n_pts = data.
da.size();
329 local_disp.setZero();
333 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
335 for (
int d = 0; d <
size(); ++d)
337 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
342 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
347 for (
long p = 0; p < n_pts; ++p)
356 const Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
357 const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
360 def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(
size(),
size());
367 Eigen::Matrix<double, n_basis, dim> gradient = delF_delU * gradient_temp.transpose();
369 G.noalias() += gradient * data.
da(p);
372 Eigen::Matrix<double, dim, n_basis> G_T = G.transpose();
374 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
375 Eigen::Matrix<double, N, 1> temp(Eigen::Map<Eigen::Matrix<double, N, 1>>(G_T.data(), G_T.size()));
382 assert(data.
x.cols() == 1);
384 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
385 const int n_pts = data.
da.size();
388 local_disp.setZero();
392 for (
size_t ii = 0; ii < bs.global.size(); ++ii)
394 for (
int d = 0; d <
size(); ++d)
396 local_disp(i, d) += bs.global[ii].val * data.
x(bs.global[ii].index *
size() + d);
401 Eigen::Matrix<double, dim, dim> def_grad(
size(),
size());
403 for (
long p = 0; p < n_pts; ++p)
412 Eigen::Matrix<double, dim, dim> jac_it = data.
vals.
jac_it[p];
413 const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
416 def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(
size(),
size());
423 Eigen::Matrix<double, dim * dim, N> delF_delU_tensor = Eigen::Matrix<double, dim * dim, N>::Zero(jac_it.size(), grad.size());
425 for (
size_t i = 0; i < local_disp.rows(); ++i)
427 for (
size_t j = 0; j < dim; ++j)
428 for (
size_t k = 0; k < dim; ++k)
429 delF_delU_tensor(dim * k + j, i * dim + j) = delF_delU(i, k);
432 Eigen::Matrix<double, N, N> hessian = delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor;
434 H += hessian * data.
da(p);
495 Eigen::MatrixXd &dstress_dmu,
496 Eigen::MatrixXd &dstress_dlambda)
const
498 const double t = data.
t;
499 const int el_id = data.
el_id;
500 const Eigen::MatrixXd &local_pts = data.
local_pts;
501 const Eigen::MatrixXd &global_pts = data.
global_pts;
502 const Eigen::MatrixXd &grad_u_i = data.
grad_u_i;
504 const Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(
size(),
size()) + grad_u_i;
507 Eigen::VectorXd sigmas;
523 Eigen::MatrixXd delJ_delF(
size(),
size());
528 delJ_delF(0, 0) = def_grad(1, 1);
529 delJ_delF(0, 1) = -def_grad(1, 0);
530 delJ_delF(1, 0) = -def_grad(0, 1);
531 delJ_delF(1, 1) = def_grad(0, 0);
533 else if (
size() == 3)
535 const Eigen::Matrix<double, 3, 1> u = def_grad.col(0);
536 const Eigen::Matrix<double, 3, 1> v = def_grad.col(1);
537 const Eigen::Matrix<double, 3, 1> w = def_grad.col(2);
539 delJ_delF.col(0) = cross<3>(v, w);
540 delJ_delF.col(1) = cross<3>(w, u);
541 delJ_delF.col(2) = cross<3>(u, v);
544 dstress_dmu = 2 * (def_grad - UVT);
545 dstress_dlambda = (sigmas.prod() - 1) * delJ_delF;
630 const double sigmaProd = sigmas.prod();
631 Eigen::Matrix<double, dim, 1> sigmaProd_noI;
632 if constexpr (dim == 2) {
633 sigmaProd_noI[0] = sigmas[1];
634 sigmaProd_noI[1] = sigmas[0];
637 sigmaProd_noI[0] = sigmas[1] * sigmas[2];
638 sigmaProd_noI[1] = sigmas[2] * sigmas[0];
639 sigmaProd_noI[2] = sigmas[0] * sigmas[1];
642 double _2mu = mu * 2;
643 Eigen::Matrix<double, dim, dim> d2E_div_dsigma2;
644 d2E_div_dsigma2.setZero();
645 d2E_div_dsigma2(0, 0) = _2mu + lambda * sigmaProd_noI[0] * sigmaProd_noI[0];
646 d2E_div_dsigma2(1, 1) = _2mu + lambda * sigmaProd_noI[1] * sigmaProd_noI[1];
647 if constexpr (dim == 3) {
648 d2E_div_dsigma2(2, 2) = _2mu + lambda * sigmaProd_noI[2] * sigmaProd_noI[2];
651 if constexpr (dim == 2) {
652 d2E_div_dsigma2(0, 1) = d2E_div_dsigma2(1, 0) = lambda * ((sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[1]);
655 d2E_div_dsigma2(0, 1) = d2E_div_dsigma2(1, 0) = lambda * (sigmas[2] * (sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[1]);
656 d2E_div_dsigma2(0, 2) = d2E_div_dsigma2(2, 0) = lambda * (sigmas[1] * (sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[2]);
657 d2E_div_dsigma2(2, 1) = d2E_div_dsigma2(1, 2) = lambda * (sigmas[0] * (sigmaProd - 1.0) + sigmaProd_noI[2] * sigmaProd_noI[1]);
660 return d2E_div_dsigma2;
708 constexpr int Cdim2 = dim * (dim - 1) / 2;
709 Eigen::Matrix<double, Cdim2, 1> BLeftCoef;
712 const double tmp = lambda / 2.0 * (sigmas.prod() - 1);
713 if constexpr (dim == 2) {
714 BLeftCoef[0] = mu - tmp;
717 BLeftCoef[0] = mu - tmp * sigmas[2];
718 BLeftCoef[1] = mu - tmp * sigmas[0];
719 BLeftCoef[2] = mu - tmp * sigmas[1];
722 Eigen::Matrix2d B[Cdim2];
723 for (
int cI = 0; cI < Cdim2; cI++) {
725 int cI_post = (cI + 1) % dim;
727 double rightCoef = dE_div_dsigma[cI] + dE_div_dsigma[cI_post];
728 double sum_sigma = sigmas[cI] + sigmas[cI_post];
729 rightCoef /= 2.0 * std::max(sum_sigma, 1.0e-12);
731 const double& leftCoef = BLeftCoef[cI];
732 B[cI](0, 0) = B[cI](1, 1) = leftCoef + rightCoef;
733 B[cI](0, 1) = B[cI](1, 0) = leftCoef - rightCoef;
737 Eigen::Matrix<double, dim * dim, dim * dim> M;
739 if constexpr (dim == 2) {
740 M(0, 0) = d2E_div_dsigma2(0, 0);
741 M(0, 3) = d2E_div_dsigma2(0, 1);
742 M.block(1, 1, 2, 2) = B[0];
743 M(3, 0) = d2E_div_dsigma2(1, 0);
744 M(3, 3) = d2E_div_dsigma2(1, 1);
748 M(0, 0) = d2E_div_dsigma2(0, 0);
749 M(0, 4) = d2E_div_dsigma2(0, 1);
750 M(0, 8) = d2E_div_dsigma2(0, 2);
751 M(4, 0) = d2E_div_dsigma2(1, 0);
752 M(4, 4) = d2E_div_dsigma2(1, 1);
753 M(4, 8) = d2E_div_dsigma2(1, 2);
754 M(8, 0) = d2E_div_dsigma2(2, 0);
755 M(8, 4) = d2E_div_dsigma2(2, 1);
756 M(8, 8) = d2E_div_dsigma2(2, 2);
758 M(1, 1) = B[0](0, 0);
759 M(1, 3) = B[0](0, 1);
760 M(3, 1) = B[0](1, 0);
761 M(3, 3) = B[0](1, 1);
763 M(5, 5) = B[1](0, 0);
764 M(5, 7) = B[1](0, 1);
765 M(7, 5) = B[1](1, 0);
766 M(7, 7) = B[1](1, 1);
768 M(2, 2) = B[2](1, 1);
769 M(2, 6) = B[2](1, 0);
770 M(6, 2) = B[2](0, 1);
771 M(6, 6) = B[2](0, 0);
775 Eigen::Matrix<double, dim*dim, dim*dim> hessian;
777 const Eigen::Matrix<double, dim, dim>& U = svd.
matrixU();
778 const Eigen::Matrix<double, dim, dim>&
V = svd.
matrixV();
779 for (
int i = 0; i < dim; i++) {
781 for (
int j = 0; j < dim; j++) {
782 int ij = i + j * dim;
783 for (
int r = 0; r < dim; r++) {
785 for (
int s = 0; s < dim; s++) {
786 int rs = r + s * dim;
792 if constexpr (dim == 2) {
793 hessian(ij, rs) = M(0, 0) * U(i, 0) *
V(j, 0) * U(r, 0) *
V(s, 0) + M(0, 3) * U(i, 0) *
V(j, 0) * U(r, 1) *
V(s, 1) + M(1, 1) * U(i, 0) *
V(j, 1) * U(r, 0) *
V(s, 1) + M(1, 2) * U(i, 0) *
V(j, 1) * U(r, 1) *
V(s, 0) + M(2, 1) * U(i, 1) *
V(j, 0) * U(r, 0) *
V(s, 1) + M(2, 2) * U(i, 1) *
V(j, 0) * U(r, 1) *
V(s, 0) + M(3, 0) * U(i, 1) *
V(j, 1) * U(r, 0) *
V(s, 0) + M(3, 3) * U(i, 1) *
V(j, 1) * U(r, 1) *
V(s, 1);
796 hessian(ij, rs) = M(0, 0) * U(i, 0) *
V(j, 0) * U(r, 0) *
V(s, 0) + M(0, 4) * U(i, 0) *
V(j, 0) * U(r, 1) *
V(s, 1) + M(0, 8) * U(i, 0) *
V(j, 0) * U(r, 2) *
V(s, 2) + M(4, 0) * U(i, 1) *
V(j, 1) * U(r, 0) *
V(s, 0) + M(4, 4) * U(i, 1) *
V(j, 1) * U(r, 1) *
V(s, 1) + M(4, 8) * U(i, 1) *
V(j, 1) * U(r, 2) *
V(s, 2) + M(8, 0) * U(i, 2) *
V(j, 2) * U(r, 0) *
V(s, 0) + M(8, 4) * U(i, 2) *
V(j, 2) * U(r, 1) *
V(s, 1) + M(8, 8) * U(i, 2) *
V(j, 2) * U(r, 2) *
V(s, 2) + M(1, 1) * U(i, 0) *
V(j, 1) * U(r, 0) *
V(s, 1) + M(1, 3) * U(i, 0) *
V(j, 1) * U(r, 1) *
V(s, 0) + M(3, 1) * U(i, 1) *
V(j, 0) * U(r, 0) *
V(s, 1) + M(3, 3) * U(i, 1) *
V(j, 0) * U(r, 1) *
V(s, 0) + M(5, 5) * U(i, 1) *
V(j, 2) * U(r, 1) *
V(s, 2) + M(5, 7) * U(i, 1) *
V(j, 2) * U(r, 2) *
V(s, 1) + M(7, 5) * U(i, 2) *
V(j, 1) * U(r, 1) *
V(s, 2) + M(7, 7) * U(i, 2) *
V(j, 1) * U(r, 2) *
V(s, 1) + M(2, 2) * U(i, 0) *
V(j, 2) * U(r, 0) *
V(s, 2) + M(2, 6) * U(i, 0) *
V(j, 2) * U(r, 2) *
V(s, 0) + M(6, 2) * U(i, 2) *
V(j, 0) * U(r, 0) *
V(s, 2) + M(6, 6) * U(i, 2) *
V(j, 0) * U(r, 2) *
V(s, 0);
800 hessian(rs, ij) = hessian(ij, rs);