Integrate grad-grad operator.
920 {
922
923 const size_t nb_row_dofs = row_data.getIndices().size();
924 const size_t nb_col_dofs = col_data.getIndices().size();
925
926 if (nb_row_dofs && nb_col_dofs) {
927
928 const bool diag = (row_data.getFieldEntities()[0]->getLocalUniqueId() ==
929 col_data.getFieldEntities()[0]->getLocalUniqueId());
930
935
936
937 double vol = OpBase::getMeasure();
938
939
940 auto t_w = OpBase::getFTensor0IntegrationWeight();
941
942
943 auto t_row_diff_base = row_data.getFTensor1DiffN<
SPACE_DIM>();
944
945
946
947 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
948 auto get_D_at_pts =
949 MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, S, DL>,
951 auto t_D_at_pts = get_D_at_pts();
952
953
954 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
955
956
958
959
960 double a = t_w * vol *
betaCoeff(t_coords(0), t_coords(1), t_coords(2));
961
962
963 int rr = 0;
965
966
967 auto t_m = OpBase::template getLocMat<SPACE_DIM>(
SPACE_DIM * rr);
968
970
971
972
973 t_rowD(
l,
j,
k) = t_D_at_pts(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
974
975
976 auto t_col_diff_base = col_data.getFTensor1DiffN<
SPACE_DIM>(gg, 0);
977
978
980 for (int cc = 0; cc <= nb_cols; ++cc) {
981
982
983 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
984
985
986 ++t_col_diff_base;
987
988
989 ++t_m;
990 }
991
992
993 ++t_row_diff_base;
994 }
995
997 ++t_row_diff_base;
998
999
1000 ++t_w;
1001 ++t_D_at_pts;
1002 ++t_coords;
1003 }
1004
1005
1006 if (diag) {
1008 auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1010 auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
1014 t_m_rr(
i,
j) = t_m_cc(
j,
i);
1015 ++t_m_rr;
1016 ++t_m_cc;
1017 }
1018 }
1019 }
1020 }
1021
1023}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
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
FTensor::Index< 'i', SPACE_DIM > i
summit Index