83 {
85
86 auto umat_op_data_ptr = boost::make_shared<UMatOpData>(this);
87 physicalPtr->matOpsDataPtr->setUserDataPtr(umat_op_data_ptr);
88
89 int nb_integration_pts = OP::getGaussPts().size2();
90
91 auto get_tag = [&]() {
93 for (
const auto &tag_range_pair : *(
physicalPtr->tagVsRangePtr)) {
94 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
95 tag_range_pair.second.end()) {
96 return tag_range_pair.first;
97 }
98 }
99 }
100#ifndef NDEBUG
103 "ADOL-C tag not found " +
105 }
106#endif
108 };
109
110 const int current_tag = get_tag();
111 auto *fe_ptr =
const_cast<FEMethod *
>(this->getFEMethod());
113
114 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->resize(DIM, DIM,
false);
115 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->resize(DIM, DIM,
false);
116 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")->resize(
117 DIM * DIM, DIM * DIM, false);
118 auto mat_grad_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"grad");
119 auto mat_P_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P");
120 auto mat_P_dF_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P_dF");
121
122#ifndef NDEBUG
123 if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
125 "Missing common data for ADOL-C evaluation");
126 }
127 if (mat_grad_ptr->size2() != DIM * DIM) {
129 "Inconsistent size of gradient matrix for ADOL-C evaluation");
130 }
131 if (mat_grad_ptr->size1() != nb_integration_pts) {
133 "Inconsistent size of gradient matrix data for ADOL-C evaluation "
134 "%zu != %d",
135 mat_grad_ptr->size1(), nb_integration_pts);
136 }
137#endif
138
140 auto get_grad_at_pts =
142 *mat_grad_ptr, nb_integration_pts);
143 auto get_P_at_pts =
145 *mat_P_ptr, nb_integration_pts);
146 auto get_P_dF_at_pts =
148 DL>::size(*mat_P_dF_ptr, nb_integration_pts);
149
153
154 auto t_grad_at_pts = get_grad_at_pts();
155 auto t_P_at_pts = get_P_at_pts();
156
157 auto next = [&]() {
158 ++t_grad_at_pts;
159 ++t_P_at_pts;
160 };
161
163
164 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
165 auto t_F = getFTensor2FromPtr<DIM, DIM>(
166 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
167 t_F(
i,
J) = t_grad_at_pts(
i,
J);
170 auto t_P = getFTensor2FromPtr<DIM, DIM>(
171 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->data().data());
172 t_P_at_pts(
i,
J) = t_P(
i,
J);
173 next();
174 }
175 }
176
182
183 auto t_grad_at_pts = get_grad_at_pts();
184 auto t_P_dF_at_pts = get_P_dF_at_pts();
186
187 auto next = [&]() {
188 ++t_grad_at_pts;
189 ++t_P_dF_at_pts;
190 };
191
192 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
193 auto t_F = getFTensor2FromPtr<DIM, DIM>(
194 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
195 t_F(
i,
J) = t_grad_at_pts(
i,
J);
198 auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
199 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")
200 ->data()
201 .data());
202 t_P_dF_at_pts(
i,
J,
k,
L) = t_P_dF(
i,
J,
k,
L);
203 next();
204 }
205 }
206
210
211 auto t_grad_at_pts = get_grad_at_pts();
212 auto next = [&]() { ++t_grad_at_pts; };
213
214 for (int gg = 0; gg != nb_integration_pts; ++gg) {
215 auto t_F = getFTensor2FromPtr<DIM, DIM>(
216 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
217 t_F(
i,
J) = t_grad_at_pts(
i,
J);
220 next();
221 }
222 }
223
225}
#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.
