18using namespace boost::numeric;
23 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
24 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
25 const std::string velocity_field,
const std::string pressure_field,
26 boost::shared_ptr<CommonData> common_data,
const EntityType type) {
29 for (
auto &sit : common_data->setOfBlocksData) {
31 if (type == MBPRISM) {
32 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
33 fe_lhs_ptr->getOpPtrVector().push_back(
35 fe_lhs_ptr->getOpPtrVector().push_back(
37 fe_lhs_ptr->getOpPtrVector().push_back(
39 fe_rhs_ptr->getOpPtrVector().push_back(
41 fe_rhs_ptr->getOpPtrVector().push_back(
43 fe_rhs_ptr->getOpPtrVector().push_back(
48 velocity_field, common_data->velPtr));
49 fe_lhs_ptr->getOpPtrVector().push_back(
51 common_data->gradVelPtr));
54 velocity_field, velocity_field, common_data, sit.second));
56 velocity_field, velocity_field, common_data, sit.second));
58 velocity_field, pressure_field, common_data, sit.second));
61 velocity_field, common_data->velPtr));
62 fe_rhs_ptr->getOpPtrVector().push_back(
64 common_data->gradVelPtr));
66 pressure_field, common_data->pressPtr));
68 fe_rhs_ptr->getOpPtrVector().push_back(
71 velocity_field, common_data, sit.second));
72 fe_rhs_ptr->getOpPtrVector().push_back(
80 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
81 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
82 const std::string velocity_field,
const std::string pressure_field,
83 boost::shared_ptr<CommonData> common_data,
const EntityType type) {
86 for (
auto &sit : common_data->setOfBlocksData) {
88 if (type == MBPRISM) {
89 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
90 fe_lhs_ptr->getOpPtrVector().push_back(
92 fe_lhs_ptr->getOpPtrVector().push_back(
94 fe_lhs_ptr->getOpPtrVector().push_back(
96 fe_rhs_ptr->getOpPtrVector().push_back(
98 fe_rhs_ptr->getOpPtrVector().push_back(
100 fe_rhs_ptr->getOpPtrVector().push_back(
105 velocity_field, velocity_field, common_data, sit.second));
107 velocity_field, pressure_field, common_data, sit.second));
109 fe_rhs_ptr->getOpPtrVector().push_back(
111 common_data->gradVelPtr));
114 pressure_field, common_data->pressPtr));
115 fe_rhs_ptr->getOpPtrVector().push_back(
117 fe_rhs_ptr->getOpPtrVector().push_back(
125 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_flux_ptr,
126 const std::string velocity_field, boost::shared_ptr<CommonData> common_data,
130 for (
auto &sit : common_data->setOfBlocksData) {
132 if (type == MBPRISM) {
133 boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
134 fe_flux_ptr->getOpPtrVector().push_back(
136 fe_flux_ptr->getOpPtrVector().push_back(
138 fe_flux_ptr->getOpPtrVector().push_back(
142 velocity_field, common_data->velPtr));
143 fe_flux_ptr->getOpPtrVector().push_back(
151 boost::shared_ptr<FaceElementForcesAndSourcesCore> dragFe,
152 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
153 std::string side_fe_name,
const std::string velocity_field,
154 const std::string pressure_field,
155 boost::shared_ptr<CommonData> common_data) {
158 auto det_ptr = boost::make_shared<VectorDouble>();
159 auto jac_ptr = boost::make_shared<MatrixDouble>();
160 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
162 for (
auto &sit : common_data->setOfFacesData) {
163 sideDragFe->getOpPtrVector().push_back(
165 common_data->gradVelPtr));
167 dragFe->getOpPtrVector().push_back(
169 dragFe->getOpPtrVector().push_back(
173 pressure_field, common_data->pressPtr));
175 dragFe->getOpPtrVector().push_back(
177 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
179 velocity_field, common_data, sit.second));
185 boost::shared_ptr<PostProcFace> postProcDragPtr,
186 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
187 std::string side_fe_name,
const std::string velocity_field,
188 const std::string pressure_field,
189 boost::shared_ptr<CommonData> common_data) {
192 auto det_ptr = boost::make_shared<VectorDouble>();
193 auto jac_ptr = boost::make_shared<MatrixDouble>();
194 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
196 for (
auto &sit : common_data->setOfFacesData) {
197 sideDragFe->getOpPtrVector().push_back(
199 common_data->gradVelPtr));
200 postProcDragPtr->getOpPtrVector().push_back(
202 postProcDragPtr->getOpPtrVector().push_back(
204 postProcDragPtr->getOpPtrVector().push_back(
207 postProcDragPtr->getOpPtrVector().push_back(
209 common_data->velPtr));
210 postProcDragPtr->getOpPtrVector().push_back(
212 common_data->pressPtr));
214 postProcDragPtr->getOpPtrVector().push_back(
216 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
217 postProcDragPtr->getOpPtrVector().push_back(
219 velocity_field, postProcDragPtr->getPostProcMesh(),
220 postProcDragPtr->getMapGaussPts(), common_data, sit.second));
248 isOnDiagonal = (row_type == col_type) && (row_side == col_side);
265 const int *row_ind = &*row_data.
getIndices().begin();
266 const int *col_ind = &*col_data.
getIndices().begin();
268 Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
269 : getFEMethod()->snes_B;
272 &*locMat.data().begin(), ADD_VALUES);
274 if (!diagonalBlock || (sYmm && !isOnDiagonal)) {
275 locMat = trans(locMat);
277 &*locMat.data().begin(), ADD_VALUES);
288 const int row_nb_base_functions = row_data.
getN().size2();
296 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
299 double w = getVolume() * getGaussPts()(3, gg);
302 for (; row_bb != row_nb_dofs / 3; row_bb++) {
304 t1(
i) = w * row_diff_base_functions(
i);
307 for (
int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
310 &locMat(3 * row_bb + 1, col_bb),
311 &locMat(3 * row_bb + 2, col_bb));
313 k(
i) -= t1(
i) * base_functions;
316 ++row_diff_base_functions;
318 for (; row_bb != row_nb_base_functions; row_bb++) {
319 ++row_diff_base_functions;
331 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
333 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
334 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
338 const int row_nb_base_functions = row_data.
getN().size2();
346 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
349 double w = getVolume() * getGaussPts()(3, gg);
350 double const alpha = w * blockData.viscousCoef;
353 for (; row_bb != row_nb_dofs / 3; row_bb++) {
357 const int final_bb = isOnDiagonal ? row_bb + 1 : col_nb_dofs / 3;
361 for (; col_bb != final_bb; col_bb++) {
363 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
365 for (
int i : {0, 1, 2}) {
366 for (
int j : {0, 1, 2}) {
368 alpha * row_diff_base_functions(
j) * col_diff_base_functions(
j);
370 alpha * row_diff_base_functions(
j) * col_diff_base_functions(
i);
375 ++col_diff_base_functions;
379 ++row_diff_base_functions;
382 for (; row_bb != row_nb_base_functions; row_bb++) {
383 ++row_diff_base_functions;
388 for (
int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
390 for (; col_bb != row_bb + 1; col_bb++) {
391 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
392 auto t_off_side = get_tensor2(locMat, 3 * col_bb, 3 * row_bb);
393 t_off_side(
i,
j) = t_assemble(
j,
i);
406 auto get_tensor2 = [](
MatrixDouble &
m,
const int r,
const int c) {
408 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 1,
c + 0),
409 &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 2,
c + 0), &
m(r + 2,
c + 1),
413 const int row_nb_base_functions = row_data.
getN().size2();
417 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
418 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
424 for (
int gg = 0; gg != row_nb_gauss_pts; gg++) {
427 double w = getVolume() * getGaussPts()(3, gg);
428 double const beta = w * blockData.inertiaCoef;
431 for (; row_bb != row_nb_dofs / 3; row_bb++) {
436 const int final_bb = col_nb_dofs / 3;
439 for (; col_bb != final_bb; col_bb++) {
441 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
443 for (
int i : {0, 1, 2}) {
444 for (
int j : {0, 1, 2}) {
446 beta * col_base_functions * t_u_grad(
i,
j) * row_base_functions;
448 beta * t_u(
j) * col_diff_base_functions(
j) * row_base_functions;
453 ++col_diff_base_functions;
454 ++col_base_functions;
458 ++row_base_functions;
461 for (; row_bb != row_nb_base_functions; row_bb++) {
462 ++row_base_functions;
480 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
481 blockData.eNts.end()) {
485 nbIntegrationPts = getGaussPts().size2();
488 CHKERR iNtegrate(row_data);
490 CHKERR aSsemble(row_data);
497 const int *indices = &*data.
getIndices().data().begin();
499 const double *vals = &*locVec.data().begin();
500 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
501 : getFEMethod()->snes_f;
514 &
v(r + 0), &
v(r + 1), &
v(r + 2));
518 locVec.resize(nbRows,
false);
522 int nb_base_functions = data.
getN().size2();
527 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
534 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
536 double w = getVolume() * getGaussPts()(3, gg);
539 const double alpha = w * blockData.viscousCoef;
541 auto t_a = get_tensor1(locVec, 0);
545 for (; rr != nbRows / 3; rr++) {
548 t_a(
i) += alpha * t_u_grad(
i,
j) * t_v_grad(
j);
549 t_a(
j) += alpha * t_u_grad(
i,
j) * t_v_grad(
i);
551 t_a(
i) -= w * t_p * t_v_grad(
i);
557 for (; rr != nb_base_functions; rr++) {
574 &
v(r + 0), &
v(r + 1), &
v(r + 2));
578 locVec.resize(nbRows,
false);
585 int nb_base_functions = data.
getN().size2();
587 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
588 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
594 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
596 double w = getVolume() * getGaussPts()(3, gg);
598 const double beta = w * blockData.inertiaCoef;
600 auto t_a = get_tensor1(locVec, 0);
604 for (; rr != nbRows / 3; rr++) {
607 t_a(
i) += beta * t_v * t_u_grad(
i,
j) * t_u(
j);
613 for (; rr != nb_base_functions; rr++) {
630 locVec.resize(nbRows,
false);
636 int nb_base_functions = data.
getN().size2();
639 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
644 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
648 double w = getVolume() * getGaussPts()(3, gg);
653 for (; rr != nbRows; rr++) {
655 for (
int ii : {0, 1, 2}) {
656 t_a -= w * t_q * t_u_grad(ii, ii);
663 for (; rr != nb_base_functions; rr++) {
677 if (type != MBVERTEX)
678 PetscFunctionReturn(0);
680 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
681 blockData.eNts.end()) {
685 const int nb_gauss_pts = getGaussPts().size2();
687 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
688 auto t_w = getFTensor0IntegrationWeight();
696 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
698 double vol = getVolume();
699 t_flux(
i) += t_w * vol * t_u(
i);
706 constexpr std::array<int, 3> indices = {0, 1, 2};
719 if (type != MBVERTEX)
720 PetscFunctionReturn(0);
722 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
723 blockData.eNts.end()) {
727 const int nb_gauss_pts = getGaussPts().size2();
729 auto pressure_drag_at_gauss_pts =
730 getFTensor1FromMat<3>(*commonData->pressureDragTract);
731 auto shear_drag_at_gauss_pts =
732 getFTensor1FromMat<3>(*commonData->shearDragTract);
733 auto total_drag_at_gauss_pts =
734 getFTensor1FromMat<3>(*commonData->totalDragTract);
742 t_pressure_drag_force(
i) = 0.0;
743 t_shear_drag_force(
i) = 0.0;
744 t_total_drag_force(
i) = 0.0;
746 auto t_w = getFTensor0IntegrationWeight();
747 auto t_normal = getFTensor1NormalsAtGaussPts();
749 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
751 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
752 double w = t_w * nrm2 * 0.5;
754 t_pressure_drag_force(
i) += w * pressure_drag_at_gauss_pts(
i);
755 t_shear_drag_force(
i) += w * shear_drag_at_gauss_pts(
i);
756 t_total_drag_force(
i) += w * total_drag_at_gauss_pts(
i);
760 ++pressure_drag_at_gauss_pts;
761 ++shear_drag_at_gauss_pts;
762 ++total_drag_at_gauss_pts;
766 constexpr std::array<int, 3> indices = {0, 1, 2};
769 &t_pressure_drag_force(0), ADD_VALUES);
771 &t_shear_drag_force(0), ADD_VALUES);
773 &t_total_drag_force(0), ADD_VALUES);
783 if (type != MBVERTEX)
784 PetscFunctionReturn(0);
786 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
787 blockData.eNts.end()) {
791 CHKERR loopSideVolumes(sideFeName, *sideFe);
793 const int nb_gauss_pts = getGaussPts().size2();
796 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
798 auto t_normal = getFTensor1NormalsAtGaussPts();
805 commonData->pressureDragTract->resize(3, nb_gauss_pts,
false);
806 commonData->pressureDragTract->clear();
808 commonData->shearDragTract->resize(3, nb_gauss_pts,
false);
809 commonData->shearDragTract->clear();
811 commonData->totalDragTract->resize(3, nb_gauss_pts,
false);
812 commonData->totalDragTract->clear();
814 auto pressure_drag_at_gauss_pts =
815 getFTensor1FromMat<3>(*commonData->pressureDragTract);
816 auto shear_drag_at_gauss_pts =
817 getFTensor1FromMat<3>(*commonData->shearDragTract);
818 auto total_drag_at_gauss_pts =
819 getFTensor1FromMat<3>(*commonData->totalDragTract);
821 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
823 double nrm2 = std::sqrt(t_normal(
i) * t_normal(
i));
824 t_unit_normal(
i) = t_normal(
i) / nrm2;
826 double mu = blockData.fluidViscosity;
828 pressure_drag_at_gauss_pts(
i) = t_p * t_unit_normal(
i);
829 shear_drag_at_gauss_pts(
i) =
830 -
mu * (t_u_grad(
i,
j) + t_u_grad(
j,
i)) * t_unit_normal(
j);
831 total_drag_at_gauss_pts(
i) =
832 pressure_drag_at_gauss_pts(
i) + shear_drag_at_gauss_pts(
i);
834 ++pressure_drag_at_gauss_pts;
835 ++shear_drag_at_gauss_pts;
836 ++total_drag_at_gauss_pts;
848 if (type != MBVERTEX)
849 PetscFunctionReturn(0);
852 bzero(def_VAL, 9 *
sizeof(
double));
854 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
855 blockData.eNts.end()) {
861 Tag th_velocity_grad;
863 Tag th_pressure_drag;
866 CHKERR postProcMesh.tag_get_handle(
"VELOCITY", 3, MB_TYPE_DOUBLE, th_velocity,
867 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
868 CHKERR postProcMesh.tag_get_handle(
"PRESSURE", 1, MB_TYPE_DOUBLE, th_pressure,
869 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
870 CHKERR postProcMesh.tag_get_handle(
"VELOCITY_GRAD", 9, MB_TYPE_DOUBLE,
872 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
874 CHKERR postProcMesh.tag_get_handle(
"PRESSURE_DRAG", 3, MB_TYPE_DOUBLE,
876 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
877 CHKERR postProcMesh.tag_get_handle(
"SHEAR_DRAG", 3, MB_TYPE_DOUBLE,
879 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
880 CHKERR postProcMesh.tag_get_handle(
"TOTAL_DRAG", 3, MB_TYPE_DOUBLE,
882 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
886 const int nb_gauss_pts = commonData->pressureDragTract->size2();
889 velGradMat.resize(3, 3);
893 pressDragVec.resize(3);
895 viscDragVec.resize(3);
897 totDragVec.resize(3);
899 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
901 for (
int i : {0, 1, 2}) {
902 for (
int j : {0, 1, 2}) {
903 velGradMat(
i,
j) = (*commonData->gradVelPtr)(
i * 3 +
j, gg);
905 velVec(
i) = (*commonData->velPtr)(
i, gg);
906 pressDragVec(
i) = (*commonData->pressureDragTract)(
i, gg);
907 viscDragVec(
i) = (*commonData->shearDragTract)(
i, gg);
908 totDragVec(
i) = (*commonData->totalDragTract)(
i, gg);
910 CHKERR postProcMesh.tag_set_data(th_velocity, &mapGaussPts[gg], 1,
912 CHKERR postProcMesh.tag_set_data(th_pressure, &mapGaussPts[gg], 1, &t_p);
913 CHKERR postProcMesh.tag_set_data(th_velocity_grad, &mapGaussPts[gg], 1,
916 CHKERR postProcMesh.tag_set_data(th_pressure_drag, &mapGaussPts[gg], 1,
919 CHKERR postProcMesh.tag_set_data(th_shear_drag, &mapGaussPts[gg], 1,
922 CHKERR postProcMesh.tag_set_data(th_total_drag, &mapGaussPts[gg], 1,
935 if (type != MBVERTEX)
936 PetscFunctionReturn(0);
938 bzero(def_VAL, 9 *
sizeof(
double));
940 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
941 blockData.eNts.end()) {
949 CHKERR postProcMesh.tag_get_handle(
"VORTICITY", 3, MB_TYPE_DOUBLE,
950 th_vorticity, MB_TAG_CREAT | MB_TAG_SPARSE,
952 CHKERR postProcMesh.tag_get_handle(
"Q", 1, MB_TYPE_DOUBLE, th_q,
953 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
954 CHKERR postProcMesh.tag_get_handle(
"L2", 1, MB_TYPE_DOUBLE, th_l2,
955 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
957 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
960 const int nb_gauss_pts = commonData->gradVelPtr->size2();
982 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
984 vorticity(0) = t_u_grad(2, 1) - t_u_grad(1, 2);
985 vorticity(1) = t_u_grad(0, 2) - t_u_grad(2, 0);
986 vorticity(2) = t_u_grad(1, 0) - t_u_grad(0, 1);
989 for (
int i = 0;
i != 3;
i++) {
990 for (
int j = 0;
j != 3;
j++) {
991 q += -0.5 * t_u_grad(
i,
j) * t_u_grad(
j,
i);
994 for (
int i = 0;
i != 3;
i++) {
995 for (
int j = 0;
j != 3;
j++) {
996 S(
i,
j) = 0.5 * (t_u_grad(
i,
j) + t_u_grad(
j,
i));
997 Omega(
i,
j) = 0.5 * (t_u_grad(
i,
j) - t_u_grad(
j,
i));
1002 for (
int i = 0;
i != 3;
i++) {
1003 for (
int j = 0;
j != 3;
j++) {
1004 for (
int k = 0;
k != 3;
k++) {
1005 M(
i,
j) += S(
i,
k) * S(
k,
j) + Omega(
i,
k) * Omega(
k,
j);
1015 int n = 3, lda = 3, info, lwork = -1;
1017 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
1018 &(eigen_values.data()[0]), &wkopt, lwork);
1020 SETERRQ1(PETSC_COMM_SELF, 1,
1021 "is something wrong with lapack_dsyev info = %d", info);
1024 info =
lapack_dsyev(
'V',
'U',
n, &(eigen_vectors.data()[0]), lda,
1025 &(eigen_values.data()[0]), work, lwork);
1027 SETERRQ1(PETSC_COMM_SELF, 1,
1028 "is something wrong with lapack_dsyev info = %d", info);
1030 map<double, int> eigen_sort;
1031 eigen_sort[eigen_values[0]] = 0;
1032 eigen_sort[eigen_values[1]] = 1;
1033 eigen_sort[eigen_values[2]] = 2;
1037 prin_vals_vect.clear();
1040 for (map<double, int>::reverse_iterator mit = eigen_sort.rbegin();
1041 mit != eigen_sort.rend(); mit++) {
1042 prin_vals_vect[ii] = eigen_values[mit->second];
1050 l2 = prin_vals_vect[1];
1061 CHKERR postProcMesh.tag_set_data(th_vorticity, &mapGaussPts[gg], 1,
1063 CHKERR postProcMesh.tag_set_data(th_q, &mapGaussPts[gg], 1, &q);
1064 CHKERR postProcMesh.tag_set_data(th_l2, &mapGaussPts[gg], 1, &l2);
1072 double r = std::sqrt(x * x + y * y + z * z);
1073 double theta = acos(x / r);
1074 double phi = atan2(y, z);
1075 double ur = cos(theta) * (1.0 + 0.5 / (r * r * r) - 1.5 / r);
1076 double ut = -sin(theta) * (1.0 - 0.25 / (r * r * r) - 0.75 / r);
1078 res[0] = ur * cos(theta) - ut * sin(theta);
1079 res[1] = ur * sin(theta) * sin(
phi) + ut * cos(theta) * sin(
phi);
1080 res[2] = ur * sin(theta) * cos(
phi) + ut * cos(theta) * cos(
phi);
VectorDouble3 stokes_flow_velocity(double x, double y, double z)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
implementation of Data Operators for Forces and Sources
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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 dofs on entity.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Calculate inverse of jacobian for face element.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Operator for fat prism element updating integration weights in the volume.
Set inverse jacobian to base functions.
Transform local reference derivatives of shape functions to global derivatives.
Assemble linear (symmetric) part of the diagonal block of the LHS Operator for assembling linear (sym...
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Assemble non-linear (non-symmetric) part of the diagonal block of the LHS Operator for assembling non...
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Assemble off-diagonal block of the LHS Operator for assembling off-diagonal block of the LHS.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode aSsemble(EntData &data)
Assemble the pressure component of the RHS vector.
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local constrains vector.
Assemble linear part of the velocity component of the RHS vector.
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local entity vector.
Assemble non-linear part of the velocity component of the RHS vector.
MoFEMErrorCode iNtegrate(EntData &data)
Calculate drag force on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
calculating volumetric flux
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Post processing output of drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFace > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
static MoFEMErrorCode setCalcVolumeFluxOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_flux_ptr, const std::string velocity_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for calculation of volume flux.