771 {
773
774 const double vol = getMeasure();
775 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*
dotUPtr);
776 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*
uPtr);
777 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*
gradUPtr);
778 auto t_h = getFTensor0FromVec(*
hPtr);
779 auto t_g = getFTensor0FromVec(*
gPtr);
780 auto t_coords = getFTensor1CoordsAtGaussPts();
781
782 auto t_row_base = row_data.getFTensor0N();
783 auto t_row_diff_base = row_data.getFTensor1DiffN<
SPACE_DIM>();
784
785 auto t_w = getFTensor0IntegrationWeight();
786
794
795 t_buoyancy_dh(
i) = 0;
797
798 for (int gg = 0; gg != nbIntegrationPts; gg++) {
799
800 const double r = t_coords(0);
802
805
806 auto t_D_dh =
get_D(2 * mu_dh);
807
808 t_inertia_force_dh(
i) = (alpha * rho_dh) * t_dot_u(
i);
809
810 t_gravity_dh(
SPACE_DIM - 1) = (alpha * rho_dh *
a0);
811 t_convection_dh(
i) = (rho_dh * alpha) * (t_u(
j) * t_grad_u(
i,
j));
812 const double t_phase_force_g_dh = -alpha *
kappa * t_g;
813 t_forces_dh(
i) = t_inertia_force_dh(
i) + t_buoyancy_dh(
i) +
814 t_gravity_dh(
i) + t_convection_dh(
i);
815
816 t_stress_dh(
i,
j) = alpha * (t_D_dh(
i,
j,
k,
l) * t_grad_u(
k,
l));
817
818 int rr = 0;
820
821 auto t_mat =
822 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr *
U_FIELD_DIM);
823 auto t_col_base = col_data.getFTensor0N(gg, 0);
824 auto t_col_diff_base = col_data.getFTensor1DiffN<
SPACE_DIM>(gg, 0);
825
826 for (int cc = 0; cc != nbCols; ++cc) {
827
828 const double bb = t_row_base * t_col_base;
829 t_mat(
i) += t_forces_dh(
i) * bb;
830 t_mat(
i) += (t_phase_force_g_dh * t_row_base) * t_col_diff_base(
i);
831 t_mat(
i) += (t_row_diff_base(
j) * t_col_base) * t_stress_dh(
i,
j);
832
834 t_mat(0) += (bb * (alpha / t_coords(0))) * (2 * mu_dh * t_u(0));
835 }
836
837 ++t_mat;
838 ++t_col_base;
839 ++t_col_diff_base;
840 }
841
842 ++t_row_base;
843 ++t_row_diff_base;
844 }
845
846 for (; rr < nbRowBaseFunctions; ++rr) {
847 ++t_row_diff_base;
848 ++t_row_base;
849 }
850
851 ++t_dot_u;
852 ++t_u;
853 ++t_grad_u;
854 ++t_h;
855 ++t_g;
856 ++t_coords;
857 ++t_w;
858 }
859
861 }
#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()
auto cylindrical
[cylindrical]
constexpr int U_FIELD_DIM
auto d_phase_function_h
Derivative of phase function with respect to h.
auto get_D
Create deviatoric stress tensor.
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
