628 {
630
631 const double vol = getMeasure();
632 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*
dotUPtr);
633 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*
uPtr);
634 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*
gradUPtr);
635 auto t_h = getFTensor0FromVec(*
hPtr);
636 auto t_g = getFTensor0FromVec(*
gPtr);
637 auto t_coords = getFTensor1CoordsAtGaussPts();
638
639 auto t_row_base = row_data.getFTensor0N();
640 auto t_row_diff_base = row_data.getFTensor1DiffN<
SPACE_DIM>();
641
642 auto t_w = getFTensor0IntegrationWeight();
643
651
652 t_buoyancy_dh(
i) = 0;
654
655 for (int gg = 0; gg != nbIntegrationPts; gg++) {
656
657 const double r = t_coords(0);
659
662
663 auto t_D_dh =
get_D(2 * mu_dh);
664
665 t_inertia_force_dh(
i) = (alpha * rho_dh) * t_dot_u(
i);
666
667 t_gravity_dh(
SPACE_DIM - 1) = (alpha * rho_dh *
a0);
668 t_convection_dh(
i) = (rho_dh * alpha) * (t_u(
j) * t_grad_u(
i,
j));
669 const double t_phase_force_g_dh = -alpha *
kappa * t_g;
670 t_forces_dh(
i) = t_inertia_force_dh(
i) + t_buoyancy_dh(
i) +
671 t_gravity_dh(
i) + t_convection_dh(
i);
672
673 t_stress_dh(
i,
j) = alpha * (t_D_dh(
i,
j,
k,
l) * t_grad_u(
k,
l));
674
675 int rr = 0;
677
678 auto t_mat =
679 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr *
U_FIELD_DIM);
680 auto t_col_base = col_data.getFTensor0N(gg, 0);
681 auto t_col_diff_base = col_data.getFTensor1DiffN<
SPACE_DIM>(gg, 0);
682
683 for (int cc = 0; cc != nbCols; ++cc) {
684
685 const double bb = t_row_base * t_col_base;
686 t_mat(
i) += t_forces_dh(
i) * bb;
687 t_mat(
i) += (t_phase_force_g_dh * t_row_base) * t_col_diff_base(
i);
688 t_mat(
i) += (t_row_diff_base(
j) * t_col_base) * t_stress_dh(
i,
j);
689
691 t_mat(0) += (bb * (alpha / t_coords(0))) * (2 * mu_dh * t_u(0));
692 }
693
694 ++t_mat;
695 ++t_col_base;
696 ++t_col_diff_base;
697 }
698
699 ++t_row_base;
700 ++t_row_diff_base;
701 }
702
703 for (; rr < nbRowBaseFunctions; ++rr) {
704 ++t_row_diff_base;
705 ++t_row_base;
706 }
707
708 ++t_dot_u;
709 ++t_u;
710 ++t_grad_u;
711 ++t_h;
712 ++t_g;
713 ++t_coords;
714 ++t_w;
715 }
716
718 }
#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 int U_FIELD_DIM
FTensor::Index< 'i', SPACE_DIM > i
constexpr CoordinateTypes coord_type
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k