52 {
54
55 int nb_integration_pts = OP::getGaussPts().size2();
56
57 constexpr int adolc_return_value = 0;
58
60
61 auto get_tag = [&]() {
62 if (physics_ptr->tagVsRangePtr) {
63 for (const auto &tag_range_pair : *physics_ptr->tagVsRangePtr) {
64 if (tag_range_pair.second.find(DomainEleOp::getFEEntityHandle()) !=
65 tag_range_pair.second.end()) {
66 return tag_range_pair.first;
67 }
68 }
69 }
70 return physics_ptr->tAg;
71 };
72
73 auto mat_grad_ptr = physics_ptr->adolcDataPtr->getCommonDataPtr("grad");
74 auto mat_P_ptr = physics_ptr->adolcDataPtr->getCommonDataPtr("P");
75 auto mat_P_dF_ptr = physics_ptr->adolcDataPtr->getCommonDataPtr("P_dF");
76
77#ifndef NDEBUG
78 if (!mat_grad_ptr || !mat_P_ptr || !mat_P_dF_ptr) {
80 "Missing common data for ADOL-C evaluation");
81 }
82 if (mat_grad_ptr->size1() != DIM * DIM) {
84 "Inconsistent size of gradient matrix for ADOL-C evaluation");
85 }
86 if (mat_grad_ptr->size2() != nb_integration_pts) {
88 "Inconsistent size of gradient matrix data for ADOL-C evaluation "
89 "%lu != %lu",
90 mat_grad_ptr->size2(), nb_integration_pts);
91 }
92#endif
93
97
98 auto t_grad_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_grad_ptr);
99 mat_P_ptr->resize(DIM * DIM, nb_integration_pts, false);
100 auto t_P_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_P_ptr);
101
102 auto next = [&]() {
103 ++t_grad_at_pts;
104 ++t_P_at_pts;
105 };
106
108
109 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
110 auto t_F = getFTensor2FromPtr<DIM, DIM>(
111 physics_ptr->adolcDataPtr->getActiveDataPtr("F")->data().data());
112 t_F(
i,
J) = t_grad_at_pts(
i,
J);
113 CHKERR physics_ptr->setActiveContinuousVector();
114 CHKERR physics_ptr->setDependentContinuousVector();
115 CHKERR physics_ptr->setParams(this->getFEEntityHandle(), gg);
116 int r = ::function(get_tag(), physics_ptr->dependentVariables.size(),
117 physics_ptr->activeVariables.size(),
118 physics_ptr->activeVariables.data(),
119 physics_ptr->dependentVariables.data());
120 if (PetscUnlikely(r < adolc_return_value)) {
122 "ADOL-C function evaluation with error");
123 }
124 CHKERR physics_ptr->getDependentContinuousVector();
125 auto t_P = getFTensor2FromPtr<DIM, DIM>(
126 physics_ptr->adolcDataPtr->getDependentDataPtr("P")->data().data());
127 t_P_at_pts(
i,
J) = t_P(
i,
J);
128 next();
129 }
130 }
131
137
138 const int number_of_active_variables = physics_ptr->activeVariables.size();
139 const int number_of_dependent_variables =
140 physics_ptr->dependentVariables.size();
141 CHKERR physics_ptr->setDependentDerivativesContinuousVector();
142 std::vector<double *> jac_ptr(number_of_dependent_variables);
143 for (
unsigned int n = 0;
n != number_of_dependent_variables; ++
n) {
145 physics_ptr
146 ->dependentVariablesDerivatives[
n * number_of_active_variables]);
147 }
148 auto t_grad_at_pts = getFTensor2FromMat<DIM, DIM>(*mat_grad_ptr);
149 mat_P_dF_ptr->resize(DIM * DIM * DIM * DIM, nb_integration_pts, false);
150 auto t_P_dF_at_pts =
151 getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(*mat_P_dF_ptr);
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 physics_ptr->adolcDataPtr->getActiveDataPtr("F")->data().data());
162 t_F(
i,
J) = t_grad_at_pts(
i,
J);
163 CHKERR physics_ptr->setActiveContinuousVector();
164 CHKERR physics_ptr->setDependentContinuousVector();
165 CHKERR physics_ptr->setParams(this->getFEEntityHandle(), gg);
166 int r = ::jacobian(get_tag(), number_of_dependent_variables,
167 number_of_active_variables,
168 physics_ptr->activeVariables.data(), jac_ptr.data());
169 if (PetscUnlikely(r < adolc_return_value)) {
171 "ADOL-C function evaluation with error");
172 }
173 CHKERR physics_ptr->getDependentDerivativesContinuousVector();
174 auto t_P_dF = getFTensor4FromPtr<DIM, DIM, DIM, DIM>(
175 physics_ptr->adolcDataPtr->getDependentDerivativesDataPtr("P_dF")
176 ->data()
177 .data());
178 t_P_dF_at_pts(
i,
J,
k,
L) = t_P_dF(
i,
J,
k,
L);
179 next();
180 }
181 }
182
184}
#define FTENSOR_INDEX(DIM, I)
#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
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'k', 3 > k
