65 {
67
68 int nb_integration_pts = OP::getGaussPts().size2();
69
70 auto get_tag = [&]() {
72 for (
const auto &tag_range_pair : *(
physicalPtr->tagVsRangePtr)) {
73 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
74 tag_range_pair.second.end()) {
75 return tag_range_pair.first;
76 }
77 }
78 }
79#ifndef NDEBUG
82 "ADOL-C tag not found " +
84 }
85#endif
87 };
88
89 const int current_tag = get_tag();
90 auto *fe_ptr =
const_cast<FEMethod *
>(this->getFEMethod());
92
93 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->resize(DIM, DIM,
false);
94 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->resize(DIM, DIM,
false);
95 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")->resize(
96 DIM * DIM, DIM * DIM, false);
97 auto mat_grad_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"grad");
98 auto mat_P_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P");
99 auto mat_P_dF_ptr =
physicalPtr->matOpsDataPtr->getCommonDataPtr(
"P_dF");
100
101#ifndef NDEBUG
102 if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
104 "Missing common data for ADOL-C evaluation");
105 }
106 if (mat_grad_ptr->size2() != DIM * DIM) {
108 "Inconsistent size of gradient matrix for ADOL-C evaluation");
109 }
110 if (mat_grad_ptr->size1() != nb_integration_pts) {
112 "Inconsistent size of gradient matrix data for ADOL-C evaluation "
113 "%zu != %d",
114 mat_grad_ptr->size1(), nb_integration_pts);
115 }
116#endif
117
119 auto get_grad_at_pts =
121 *mat_grad_ptr, nb_integration_pts);
122 auto get_P_at_pts =
124 *mat_P_ptr, nb_integration_pts);
125 auto get_P_dF_at_pts =
127 DL>
::size(*mat_P_dF_ptr, nb_integration_pts);
128
132
133 auto t_grad_at_pts = get_grad_at_pts();
134 auto t_P_at_pts = get_P_at_pts();
135
136 auto next = [&]() {
137 ++t_grad_at_pts;
138 ++t_P_at_pts;
139 };
140
141 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
142 auto t_F = getFTensor2FromPtr<DIM, DIM>(
143 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
144 t_F(
i,
J) = t_grad_at_pts(
i,
J);
147 auto t_P = getFTensor2FromPtr<DIM, DIM>(
148 physicalPtr->matOpsDataPtr->getDependentDataPtr(
"P")->data().data());
149 t_P_at_pts(
i,
J) = t_P(
i,
J);
150 next();
151 }
152 }
153
159
160 auto t_grad_at_pts = get_grad_at_pts();
161 auto t_P_dF_at_pts = get_P_dF_at_pts();
162 auto next = [&]() {
163 ++t_grad_at_pts;
164 ++t_P_dF_at_pts;
165 };
166
167 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
168 auto t_F = getFTensor2FromPtr<DIM, DIM>(
169 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
170 t_F(
i,
J) = t_grad_at_pts(
i,
J);
173 auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
174 physicalPtr->matOpsDataPtr->getDependentDerivativesDataPtr(
"P_dF")
175 ->data()
176 .data());
177 t_P_dF_at_pts(
i,
J,
k,
L) = t_P_dF(
i,
J,
k,
L);
178 next();
179 }
180 }
181
185
186 auto t_grad_at_pts = get_grad_at_pts();
187 auto next = [&]() { ++t_grad_at_pts; };
188
189 for (int gg = 0; gg != nb_integration_pts; ++gg) {
190 auto t_F = getFTensor2FromPtr<DIM, DIM>(
191 physicalPtr->matOpsDataPtr->getActiveDataPtr(
"F")->data().data());
192 t_F(
i,
J) = t_grad_at_pts(
i,
J);
195 next();
196 }
197 }
198
200}
#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.
