205 {
207
213
214 int nb_integration_pts = row_data.
getN().size1();
215 int row_nb_dofs = row_data.
getIndices().size();
216 int col_nb_dofs = col_data.
getIndices().size();
217
221
222 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
223 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
224
225 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
226 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
227
228 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
229 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
230
231 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
232 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
233
234 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
235 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
236
237 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
238 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
239
240 );
241 };
242
249
251
252 auto get_dP = [&]() {
253 auto get_tensor_for_dP =
255 DL>::size(dP, nb_integration_pts);
256
257 auto mat_ops_data =
physicsPtr->matOpsDataPtr;
258
260 mat_ops_data->getCommonDataPtr("PAtPts"));
261 auto t_P_du = getFTensor4FromMat<3, 3, 3, 3, -1,
DL>(
262 mat_ops_data->getCommonDataPtr("PAtPts_du"));
264 mat_ops_data->getCommonDataPtr("adjointPdstretchAtPts"));
266 mat_ops_data->getCommonDataPtr("stretchH1AtPts"));
267 auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1,
DL>(
268 mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
270 mat_ops_data->getCommonDataPtr("wGradH1AtPts"));
272 mat_ops_data->getCommonDataPtr("eigenVals"));
274 mat_ops_data->getCommonDataPtr("eigenVecs"));
275
278
279 auto t_dP = get_tensor_for_dP();
280 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
281
283
288 break;
289 default:
291 }
292
294 t_Ldiff_u(
i,
j, L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n, L);
296 -t_Ldiff_u(
i,
j, L) * (t_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
297 t_dP(L,
J) += (ts_a *
alphaU) *
298 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
299
304 t_approx_P_adjoint_dstretch(
i,
j) - t_P(
i,
k) * t_h1(
j,
k);
306 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
307 t_deltaP_sym(
i,
j) /= 2.0;
308 auto nb_uniq =
314 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
315 }
316
317 ++t_P;
318 ++t_dP;
319 ++t_P_du;
320 ++t_approx_P_adjoint_dstretch;
321 ++t_u;
322 ++t_diff_u;
323 ++t_grad_h1;
324 ++t_eigen_vals;
325 ++t_eigen_vecs;
326 }
327
328 return get_tensor_for_dP();
329 };
330
331 int row_nb_base_functions = row_data.
getN().size2();
334
335 auto t_dP = get_dP();
339
340 for (int gg = 0; gg != nb_integration_pts; ++gg) {
342
343 int rr = 0;
344 for (; rr != row_nb_dofs /
size_symm; ++rr) {
347 auto t_m = get_ftensor2(
K,
size_symm * rr, 0);
348 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
349 const double b =
a * t_row_base_fun * t_col_base_fun;
350 t_m(L,
J) -= b * t_dP(L,
J);
352 (t_row_grad_fun(
i) * t_col_grad_fun(
i));
353 t_m(L,
J) +=
c * t_kd_sym(L,
J);
354 ++t_m;
355 ++t_col_base_fun;
356 ++t_col_grad_fun;
357 }
358 ++t_row_base_fun;
359 ++t_row_grad_fun;
360 }
361
362 for (; rr != row_nb_base_functions; ++rr) {
363 ++t_row_base_fun;
364 ++t_row_grad_fun;
365 }
366
367 ++t_w;
368 ++t_dP;
369 }
371}
#define FTENSOR_INDEX(DIM, I)
Fourth-order differential tensor symmetric in both index pairs.
Kronecker Delta class symmetric.
Mapping from symmetric tensor indices to packed storage index.
#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 c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'J', DIM1 > J
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
static constexpr auto size_symm
UBlasMatrix< double > MatrixDouble
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
FTensor::Index< 'm', 3 > m
static enum StretchSelector stretchSelector
static enum RotSelector gradApproximator
static boost::function< double(const double)> f
static boost::function< double(const double)> dd_f
static boost::function< double(const double)> d_f
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getVolume() const
element volume (linear geometry)
MatrixDouble K
local tangent matrix
Public Member Functions inherited from OpAssembleVolume