704 {
706
707
709 const auto in_the_loop =
711
712 auto not_side = [](auto s) {
714 };
715
718 &*base_mat.data().begin());
719 };
720
721 if (in_the_loop > 0) {
722
723
724 auto t_normal = getFTensor1Normal();
725 const auto nb_gauss_pts = getGaussPts().size2();
726
728
729
730 const auto nb_rows =
sideDataPtr->indicesRowSideMap[s0].size();
731
732 if (nb_rows) {
733
736
737
738 const auto opposite_s0 = not_side(s0);
740#ifndef NDEBUG
741 const auto opposite_sense_row =
sideDataPtr->senseMap[opposite_s0];
742 if (sense_row * opposite_sense_row > 0)
744 "Should be opposite sign");
745#endif
746
747
748 const auto nb_row_base_functions =
750
751 auto t_w = getFTensor0IntegrationWeight();
758
759 auto next = [&]() {
760 for (auto &t_l : arr_t_l)
761 ++t_l;
762 for (auto &t_vel : arr_t_vel)
763 ++t_vel;
764 };
765
766#ifndef NDEBUG
767 if (nb_gauss_pts !=
sideDataPtr->rowBaseSideMap[s0].size1())
769 "Inconsistent number of DOFs");
770#endif
771
773 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
776 const auto dot = sense_row * (t_normal(
i) * t_vel(
i));
777 const auto l_upwind_side = (dot > 0) ? s0 : opposite_s0;
778 const auto l_upwind = arr_t_l[l_upwind_side];
780 t_res(
I,
J) = t_w * dot * l_upwind(
I,
J);
781 next();
782 ++t_w;
783
784 auto t_res_skeleton =
786 auto rr = 0;
787 for (; rr != nb_rows; ++rr) {
788 t_res_skeleton(
I,
J) += t_row_base * t_res(
I,
J);
789 ++t_row_base;
790 ++t_res_skeleton;
791 }
792 for (; rr < nb_row_base_functions; ++rr) {
793 ++t_row_base;
794 }
795 }
796
797 CHKERR ::VecSetValues(getTSf(),
801 }
802 }
803 }
804
806}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'I', DIM1 > I
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
constexpr auto make_array(Arg &&...arg)
Create Array.
auto get_ntensor(T &base_mat)