944 {
946
949
950 auto t_normal = getFTensor1Normal();
951 const double nrm2 = sqrt(t_normal(
i) * t_normal(
i));
953 t_unit_normal(
i) = t_normal(
i) / nrm2;
955 int nb_integration_pts = data.
getN().size1();
956 int nb_base_functions = data.
getN().size2() / 3;
957 double ts_t = getFEMethod()->ts_t;
958
959 auto integrate_matrix = [&]() {
961
962 auto t_w = getFTensor0IntegrationWeight();
963 matPP.resize(nb_dofs / 3, nb_dofs / 3,
false);
965
967 for (int gg = 0; gg != nb_integration_pts; ++gg) {
968
969 int rr = 0;
970 for (; rr != nb_dofs / 3; ++rr) {
971 const double a = t_row_base_fun(
i) * t_unit_normal(
i);
973 for (int cc = 0; cc != nb_dofs / 3; ++cc) {
974 const double b = t_col_base_fun(
i) * t_unit_normal(
i);
976 ++t_col_base_fun;
977 }
978 ++t_row_base_fun;
979 }
980
981 for (; rr != nb_base_functions; ++rr)
982 ++t_row_base_fun;
983
984 ++t_w;
985 }
986
988 };
989
990 auto integrate_rhs = [&](auto &bc) {
992
993 auto t_w = getFTensor0IntegrationWeight();
994 vecPv.resize(3, nb_dofs / 3,
false);
996
998 double ts_t = getFEMethod()->ts_t;
999
1000 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1001 int rr = 0;
1002 for (; rr != nb_dofs / 3; ++rr) {
1003 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1004 for (
int dd = 0;
dd != 3; ++
dd)
1005 if (bc.flags[dd])
1007 ++t_row_base_fun;
1008 }
1009 for (; rr != nb_base_functions; ++rr)
1010 ++t_row_base_fun;
1011 ++t_w;
1012 }
1014 };
1015
1016 auto integrate_rhs_cook = [&](auto &bc) {
1018
1019 vecPv.resize(3, nb_dofs / 3,
false);
1021
1022 auto t_w = getFTensor0IntegrationWeight();
1024 auto t_coords = getFTensor1CoordsAtGaussPts();
1025
1026 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1027
1028 auto calc_tau = [](double y) {
1029 y -= 44;
1030 y /= (60 - 44);
1031 return -y * (y - 1) / 0.25;
1032 };
1033
1034 const double tau = calc_tau(t_coords(1));
1035
1036 int rr = 0;
1037 for (; rr != nb_dofs / 3; ++rr) {
1038 const double t = ts_t * t_w * t_row_base_fun(
i) * t_unit_normal(
i);
1039 for (
int dd = 0;
dd != 3; ++
dd)
1040 if (bc.flags[dd])
1041 vecPv(dd, rr) +=
t * tau * bc.vals[
dd];
1042 ++t_row_base_fun;
1043 }
1044
1045 for (; rr != nb_base_functions; ++rr)
1046 ++t_row_base_fun;
1047 ++t_w;
1048 ++t_coords;
1049 }
1051 };
1052
1053
1056 if (bc.faces.find(fe_ent) != bc.faces.end()) {
1057
1059 if (nb_dofs) {
1060
1061 CHKERR integrate_matrix();
1062 if (std::regex_match(bc.blockName, std::regex(".*COOK.*")))
1063 CHKERR integrate_rhs_cook(bc);
1064 else
1065 CHKERR integrate_rhs(bc);
1066
1067 const auto info =
1069 nb_dofs / 3, &*
vecPv.data().begin(), nb_dofs / 3);
1070 if (info > 0)
1072 "The leading minor of order %i is "
1073 "not positive; definite;\nthe "
1074 "solution could not be computed",
1075 info);
1076
1077 for (
int dd = 0;
dd != 3; ++
dd)
1078 if (bc.flags[dd])
1079 for (int rr = 0; rr != nb_dofs / 3; ++rr)
1081 }
1082 }
1083 }
1084
1086}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'i', SPACE_DIM > i
static __CLPK_integer lapack_dposv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *b, __CLPK_integer ldb)
FTensor::Index< 'j', 3 > j
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
constexpr double t
plate stiffness
boost::shared_ptr< TractionBcVec > bcData
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity