Integrate grad-grad operator.
180 {
182
183 const double vol = getMeasure();
184 auto t_w = getFTensor0IntegrationWeight();
185 auto t_row_base = row_data.getFTensor0N();
186 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
187 auto t_coords = getFTensor1CoordsAtGaussPts();
188 auto t_normal = getFTensor1NormalsAtGaussPts();
189
190 auto t_u = getFTensor1FromMat<3>(*
uPtr);
191 auto t_grad_u = getFTensor2FromMat<3, 3>(*
uGradPtr);
192
193 auto get_t_mat = [&](const int rr) {
196
198
200 };
201
202 const auto ts_a = getFEMethod()->ts_a;
203
205
206 const auto a = std::sqrt(t_coords(
i) * t_coords(
i));
207 const auto sin_fi = t_coords(2) /
a;
208 const auto f = 2 *
omega * sin_fi;
209
211 t_r(
i) = t_normal(
i);
213
216 t_P(
i,
j) = t_r(
i) * t_r(
j);
218
221
224 t_Q(
m,
i) * (ts_a *
t_kd(
i,
j) + t_grad_u(
i,
j) +
f * t_A(
i,
j)) +
226
227 const double alpha = t_w * vol;
228
229 int rr = 0;
230 for (; rr !=
nbRows / 3; rr++) {
231
232 auto t_col_base = col_data.getFTensor0N(gg, 0);
233 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
234 auto t_mat = get_t_mat(3 * rr);
235
236 for (
int cc = 0; cc !=
nbCols / 3; cc++) {
237 t_mat(
i,
j) += (alpha * t_row_base * t_col_base) * t_rhs_du(
i,
j);
238 t_mat(
i,
j) += (alpha *
mu) * t_Q(
i,
j) *
239 (t_row_diff_base(
m) * t_col_diff_base(
m));
240 ++t_col_diff_base;
241 ++t_col_base;
242 ++t_mat;
243 }
244 ++t_row_base;
245 ++t_row_diff_base;
246 }
248 ++t_row_base;
249 ++t_row_diff_base;
250 }
251
252 ++t_w;
253 ++t_coords;
254 ++t_normal;
255 ++t_u;
256 ++t_grad_u;
257 }
258
260 }
Tensor1< T, Tensor_Dim > normalize()
#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()
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'm', 3 > m
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