412 {
414
417 auto t_L = FTensor::SymmLTensor<double, 3>();
418 auto t_diff = FTensor::DiffTensor<double>();
420
421 int nb_integration_pts = row_data.
getN().size1();
422 int row_nb_dofs = row_data.
getIndices().size();
423 int col_nb_dofs = col_data.
getIndices().size();
424
428
429 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
430 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
431
432 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
433 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
434
435 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
436 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
437
438 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
439 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
440
441 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
442 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
443
444 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
445 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
446
447 );
448 };
449
456
458
459 auto get_dP = [&]() {
460 auto get_tensor_for_dP =
462 DL>::size(dP, nb_integration_pts);
463
464 auto mat_ops_data =
physicsPtr->matOpsDataPtr;
465
467 mat_ops_data->getCommonDataPtr("PAtPts"));
468 auto t_P_du = getFTensor4FromMat<3, 3, 3, 3, -1,
DL>(
469 mat_ops_data->getCommonDataPtr("PAtPts_du"));
471 mat_ops_data->getCommonDataPtr("adjointPdstretchAtPts"));
473 mat_ops_data->getCommonDataPtr("stretchH1AtPts"));
474 auto t_diff_u = getFTensor4FromMat<3, 3, 3, 3, -1,
DL>(
475 mat_ops_data->getCommonDataPtr("diffStretchH1AtPts"));
477 mat_ops_data->getCommonDataPtr("wGradH1AtPts"));
479 mat_ops_data->getCommonDataPtr("eigenVals"));
481 mat_ops_data->getCommonDataPtr("eigenVecs"));
483 mat_ops_data->getCommonDataPtr("invPlasticF"));
484
487
488 auto t_dP = get_tensor_for_dP();
489 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
490
492
496 t_h1(
i,
j) = t_grad_h1(
i,
k) * t_invPlasticF(
k,
j) +
t_kd(
i,
j);
497 break;
498 default:
500 }
501
503 t_Ldiff_u(
i,
j, L) = t_diff_u(
i,
j,
m,
n) * t_L(
m,
n, L);
505 -t_Ldiff_u(
i,
j, L) * (t_P_du(
i,
j,
k,
l) * t_Ldiff_u(
k,
l,
J));
506 t_dP(L,
J) += (ts_a *
alphaU) *
507 (t_L(
i,
j, L) * (t_diff(
i,
j,
k,
l) * t_L(
k,
l,
J)));
508
513 t_approx_P_adjoint_dstretch(
i,
j) - t_P(
i,
k) * t_h1(
j,
k);
515 t_deltaP_sym(
i,
j) = (t_deltaP(
i,
j) || t_deltaP(
j,
i));
516 t_deltaP_sym(
i,
j) /= 2.0;
522 t_dP(L,
J) += t_L(
i,
j, L) * (t_diff2_uP2(
i,
j,
k,
l) * t_L(
k,
l,
J));
523 }
524
525 ++t_P;
526 ++t_dP;
527 ++t_P_du;
528 ++t_approx_P_adjoint_dstretch;
529 ++t_u;
530 ++t_diff_u;
531 ++t_grad_h1;
532 ++t_eigen_vals;
533 ++t_eigen_vecs;
534 ++t_invPlasticF;
535 }
536
537 return get_tensor_for_dP();
538 };
539
540 int row_nb_base_functions = row_data.
getN().size2();
543
544 auto t_dP = get_dP();
548 auto mat_ops_data =
physicsPtr->matOpsDataPtr;
549 auto t_plasticF = getFTensor2SymmetricFromMat<3, -1,
DL>(
550 mat_ops_data->getCommonDataPtr("plasticF"));
552 mat_ops_data->getCommonDataPtr("invPlasticF"));
553
554 for (int gg = 0; gg != nb_integration_pts; ++gg) {
556 const double a =
v * t_w * det_plasticF;
557
558 int rr = 0;
559 for (; rr != row_nb_dofs /
size_symm; ++rr) {
561 t_row_grad_intermediate(
j) = t_row_grad_fun(
i) * t_invPlasticF(
i,
j);
564 auto t_m = get_ftensor2(
K,
size_symm * rr, 0);
565 for (
int cc = 0; cc != col_nb_dofs /
size_symm; ++cc) {
567 t_col_grad_intermediate(
j) = t_col_grad_fun(
i) * t_invPlasticF(
i,
j);
568 const double b =
a * t_row_base_fun * t_col_base_fun;
569 t_m(L,
J) -= b * t_dP(L,
J);
570 const double c = (
a *
alphaGradU * ts_a) * (t_row_grad_intermediate(
j) *
571 t_col_grad_intermediate(
j));
572 t_m(L,
J) +=
c * t_kd_sym(L,
J);
573 ++t_m;
574 ++t_col_base_fun;
575 ++t_col_grad_fun;
576 }
577 ++t_row_base_fun;
578 ++t_row_grad_fun;
579 }
580
581 for (; rr != row_nb_base_functions; ++rr) {
582 ++t_row_base_fun;
583 ++t_row_grad_fun;
584 }
585
586 ++t_w;
587 ++t_dP;
588 ++t_plasticF;
589 ++t_invPlasticF;
590 }
592}
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#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
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
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