1588 {
1590
1595
1598
1600
1602
1603 auto nb_gauss_pts = getGaussPts().size2();
1604 auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
1605 auto objective_dstress =
1606 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1607 auto objective_dstrain =
1608 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1609 auto obj_grad = boost::make_shared<MatrixDouble>(
SPACE_DIM, nb_gauss_pts);
1610
1611 auto evaluate_python = [&]() {
1613 auto &coords = OP::getCoordsAtGaussPts();
1616 objective_ptr);
1619 objective_dstress);
1622 objective_dstrain);
1623
1624 auto t_grad_u =
1626 auto t_D =
1632
1634 auto t_obj_dstress =
1636 auto t_obj_dstrain =
1638
1639 auto vol = OP::getMeasure();
1640 auto t_w = getFTensor0IntegrationWeight();
1641 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1642
1646
1648 t_diff_inv_jac(
i,
j) =
1649 -(t_inv_jac(
i,
I) * t_diff_jac(
I,
J)) * t_inv_jac(
J,
j);
1651 t_diff_grad(
i,
j) = t_grad_u(
i,
k) * t_diff_inv_jac(
k,
j);
1652
1654 t_d_strain(
i,
j) = t_diff_symm(
i,
j,
k,
l) * (
1655
1657
1658 +
1659
1661
1662 );
1663
1664 auto alpha = t_w * vol;
1665
1666 (*globObjectivePtr) += alpha * t_obj;
1667 (*globObjectiveGradPtr) +=
1668 alpha *
1669 (
1670
1671 t_obj_dstress(
i,
j) * (t_D(
i,
j,
k,
l) * t_d_strain(
k,
l))
1672
1673 +
1674
1675 t_obj_dstrain(
i,
j) * t_d_strain(
i,
j)
1676
1677 +
1678
1679 t_obj * t_cof
1680
1681 );
1682
1683 ++t_w;
1684 ++t_jac;
1685 ++t_diff_jac;
1686 ++t_cof;
1687
1688 ++t_obj;
1689 ++t_obj_dstress;
1690 ++t_obj_dstrain;
1691
1692 ++t_grad_u;
1693 ++t_d_grad;
1694 }
1696 };
1697
1698 CHKERR evaluate_python();
1699
1701 }
#define FTENSOR_INDEX(DIM, I)
auto diff_symmetrize(FTensor::Number< DIM >)
constexpr int SPACE_DIM
[Define dimension]
constexpr IntegrationType I
#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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
