1506 {
1512
1514 auto nb_gauss_pts = getGaussPts().size2();
1515
1516 auto objective_dstress =
1517 boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
1518 auto objective_dstrain =
1519 boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
1520 auto objective_du =
1521 boost::make_shared<MatrixDouble>(nb_gauss_pts,
SPACE_DIM);
1522
1523 auto evaluate_python = [&]() {
1525 auto &coords = OP::getCoordsAtGaussPts();
1528 objective_dstress);
1531 objective_dstrain);
1534 objective_du);
1535
1536 auto vol = OP::getMeasure();
1537 auto t_w = OP::getFTensor0IntegrationWeight();
1538
1539 auto t_D =
1540 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(
commPtr->matDPtr));
1543
1544 auto t_obj_dstress =
1545 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
1546 auto t_obj_dstrain =
1547 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
1548 auto t_obj_du = getFTensor1FromMat<SPACE_DIM>(*objective_du);
1549
1550 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1551 const double alpha = t_w * vol;
1553 t_adjoint_stress(
i,
j) =
1554 t_D(
i,
j,
k,
l) * t_obj_dstress(
k,
l) + t_obj_dstrain(
i,
j);
1555
1556 auto t_nf = OP::template getNf<SPACE_DIM>();
1557 int rr = 0;
1558 for (; rr != OP::nbRows /
SPACE_DIM; rr++) {
1559 t_nf(
j) += alpha * t_row_grad(
i) * t_adjoint_stress(
i,
j);
1560 t_nf(
j) += alpha * t_row_base * t_obj_du(
j);
1561
1562 ++t_row_grad;
1563 ++t_row_base;
1564 ++t_nf;
1565 }
1566
1567 for (; rr < OP::nbRowBaseFunctions; ++rr) {
1568 ++t_row_grad;
1569 ++t_row_base;
1570 }
1571 ++t_obj_dstrain;
1572 ++t_obj_dstress;
1573 ++t_obj_du;
1574 ++t_w;
1575 }
1577 };
1578 CHKERR evaluate_python();
1580 }
#define FTENSOR_INDEX(DIM, 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.
constexpr int SPACE_DIM
[Define dimension]
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.