Integrate grad-grad operator.
1145 {
1147
1149
1156
1157 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
1158 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
1159
1160#ifndef NDEBUG
1163 "Wrong number of integration pts (%zu != %d)", row_hessian.size1(),
1165 }
1166 if (row_hessian.size2() !=
1169 "Wrong number of base functions (%zu != %d)",
1172 }
1175 "Wrong number of base functions (%zu < %d)", row_hessian.size2(),
1177 }
1180 "Wrong number of integration pts (%zu != %d)", col_hessian.size1(),
1182 }
1185 "Wrong number of base functions (%zu < %d)", col_hessian.size2(),
1187 }
1188#endif
1189
1190
1191 double vol = OpBase::getMeasure();
1192
1193
1194 auto t_w = OpBase::getFTensor0IntegrationWeight();
1195
1196 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1197 &*row_hessian.data().begin());
1198
1199
1200 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
1201 auto get_D_at_pts = MatrixSizeHelper<
1205 auto t_D_at_pts = get_D_at_pts();
1206
1208
1209
1211
1212
1213 double a = t_w * vol;
1214
1215
1217
1218
1219 int rr = 0;
1221
1223 t_rowD(
k,
l) = t_D_at_pts(
i,
j,
k,
l) * (
a * t_row_diff2(
i,
j));
1224
1225
1226 auto t_col_diff2 =
1227 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1228
1229 auto t_mat = getFTensor2FromPtr<FIELD_DIM, FIELD_DIM, FIELD_DIM>(
1231
1232
1234
1235
1236 t_mat(
I,
J) += (t_rowD(
i,
j) * t_col_diff2(
i,
j)) *
t_kd(
I,
J);
1237
1238
1239 ++t_col_diff2;
1240
1241
1242 ++t_mat;
1243 }
1244
1245
1246 ++t_row_diff2;
1247 }
1248
1250 ++t_row_diff2;
1251
1252
1253 ++t_w;
1254 ++t_D_at_pts;
1255 }
1256 }
1257
1259}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > 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
constexpr IntegrationType I
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions