46 {
48
55
56 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
58
59 auto &
locMat = AssemblyDomainEleOp::locMat;
60
61 const size_t nb_integration_pts = row_data.getN().size1();
62 const size_t nb_row_base_functions = row_data.getN().size2();
63
64 if (AssemblyDomainEleOp::rowType == MBVERTEX &&
65 AssemblyDomainEleOp::rowSide == 0) {
67 auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(
resDiff);
68 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*
mDPtr);
69 auto t_logC_dC =
71 auto t_grad =
75 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
79 t_DLogC_dC(
i,
j,
k,
l) = t_D(
m,
n,
k,
l) * t_logC_dC(
m,
n,
i,
j);
81 t_FDLogC_dC(
i,
j,
k,
l) = t_F(
i,
m) * t_DLogC_dC(
m,
j,
k,
l);
83 t_res_diff(
i,
j, L) = t_FDLogC_dC(
i,
j,
k,
l) * t_L(
k,
l, L);
84 ++t_logC_dC;
85 ++t_grad;
86 ++t_res_diff;
87 }
88 }
89
90 auto t_w = AssemblyDomainEleOp::getFTensor0IntegrationWeight();
91 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
92 auto t_res_diff = getFTensor3FromMat<DIM, DIM, size_symm>(
resDiff);
93
94 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
95 double alpha = AssemblyDomainEleOp::getMeasure() * t_w;
96
97 size_t rr = 0;
98 for (; rr != AssemblyDomainEleOp::nbRows / DIM; ++rr) {
99
100 auto t_mat =
102
104 t_tmp(
i, L) = (t_res_diff(
i,
j, L) * (alpha * t_row_diff_base(
j)));
105
106 auto t_col_base = col_data.getFTensor0N(gg, 0);
107 for (
size_t cc = 0; cc != AssemblyDomainEleOp::nbCols /
size_symm; ++cc) {
108
109 t_mat(
i, L) -= (t_col_base * t_tmp(
i, L));
110
111 ++t_mat;
112 ++t_col_base;
113 }
114
115 ++t_row_diff_base;
116 }
117
118 for (; rr < nb_row_base_functions; ++rr)
119 ++t_row_diff_base;
120
121 ++t_w;
122 ++t_res_diff;
123 }
124
126}
#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< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto symm_L_tensor(FTensor::Number< DIM >)
static FTensor::Tensor2< FTensor::PackPtr< double *, 3 >, 2, 3 > get_mat_vector_dtensor_sym(size_t rr, MatrixDouble &mat, FTensor::Number< 2 >)
FTensor::Index< 'm', 3 > m