53 {
55
56 int nb_integration_pts = OP::getGaussPts().size2();
57
58 auto get_tag = [&]() {
60 for (
const auto &tag_range_pair : *(
physicalPtr->tagVsRangePtr)) {
61 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
62 tag_range_pair.second.end()) {
63 return tag_range_pair.first;
64 }
65 }
66 }
67#ifndef NDEBUG
70 "ADOL-C tag not found " +
72 }
73#endif
75 };
76
77 const int current_tag = get_tag();
78 auto *fe_ptr =
const_cast<FEMethod *
>(this->getFEMethod());
80
81 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->resize(DIM, DIM,
false);
82 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->resize(DIM, DIM,
false);
83 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")->resize(
84 DIM * DIM, DIM * DIM, false);
85 auto mat_grad_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"grad");
86 auto mat_P_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P");
87 auto mat_P_dF_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P_dF");
88
89#ifndef NDEBUG
90 if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
92 "Missing common data for ADOL-C evaluation");
93 }
94 if (mat_grad_ptr->size2() != DIM * DIM) {
96 "Inconsistent size of gradient matrix for ADOL-C evaluation");
97 }
98 if (mat_grad_ptr->size1() != nb_integration_pts) {
100 "Inconsistent size of gradient matrix data for ADOL-C evaluation "
101 "%zu != %d",
102 mat_grad_ptr->size1(), nb_integration_pts);
103 }
104#endif
105
107 auto get_grad_at_pts =
109 *mat_grad_ptr, nb_integration_pts);
110 auto get_P_at_pts =
112 *mat_P_ptr, nb_integration_pts);
113 auto get_P_dF_at_pts =
115 DL>::size(*mat_P_dF_ptr, nb_integration_pts);
116
120
121 auto t_grad_at_pts = get_grad_at_pts();
122 auto t_P_at_pts = get_P_at_pts();
123
124 auto next = [&]() {
125 ++t_grad_at_pts;
126 ++t_P_at_pts;
127 };
128
130
131 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
132 auto t_F = getFTensor2FromPtr<DIM, DIM>(
133 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
134 t_F(
i,
J) = t_grad_at_pts(
i,
J);
137 auto t_P = getFTensor2FromPtr<DIM, DIM>(
138 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->data().data());
139 t_P_at_pts(
i,
J) = t_P(
i,
J);
140 next();
141 }
142 }
143
149
150 auto t_grad_at_pts = get_grad_at_pts();
151 auto t_P_dF_at_pts = get_P_dF_at_pts();
153
154 auto next = [&]() {
155 ++t_grad_at_pts;
156 ++t_P_dF_at_pts;
157 };
158
159 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
160 auto t_F = getFTensor2FromPtr<DIM, DIM>(
161 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
162 t_F(
i,
J) = t_grad_at_pts(
i,
J);
165 auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
166 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")
167 ->data()
168 .data());
169 t_P_dF_at_pts(
i,
J,
k,
L) = t_P_dF(
i,
J,
k,
L);
170 next();
171 }
172 }
173
177
178 auto t_grad_at_pts = get_grad_at_pts();
179 auto next = [&]() { ++t_grad_at_pts; };
180
181 for (int gg = 0; gg != nb_integration_pts; ++gg) {
182 auto t_F = getFTensor2FromPtr<DIM, DIM>(
183 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
184 t_F(
i,
J) = t_grad_at_pts(
i,
J);
187 next();
188 }
189 }
190
192}
#define FTENSOR_INDEX(DIM, I)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#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< 'k', 3 > k
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
Structure for user loop methods on finite elements.
EntityHandle getFEEntityHandle() const
Get the entity handle of the current finite element.
