965 {
967
968#ifndef NDEBUG
971 "DataAtIntegrationPts pointer is null");
972#endif
973
975 if (!nb_dofs)
977
978 const int nb_integration_pts = data.
getN().size1();
979
980#ifndef NDEBUG
981 if (
dataAtPts->divPAtPts.size2() != nb_integration_pts)
983 "divPAtPts columns (%d) != nb integration points (%d)",
984 static_cast<int>(
dataAtPts->divPAtPts.size2()),
985 nb_integration_pts);
986 if (
dataAtPts->wL2DotAtPts.size2() != nb_integration_pts)
988 "wL2DotAtPts columns (%d) != nb integration points (%d)",
989 static_cast<int>(
dataAtPts->wL2DotAtPts.size2()),
990 nb_integration_pts);
991 if (
dataAtPts->varWL2.size2() != nb_integration_pts)
993 "varWL2 columns (%d) != nb integration points (%d)",
994 static_cast<int>(
dataAtPts->varWL2.size2()), nb_integration_pts);
995 if (
dataAtPts->varRotAxis.size2() != nb_integration_pts)
997 "varRotAxis columns (%d) != nb integration points (%d)",
998 static_cast<int>(
dataAtPts->varRotAxis.size2()),
999 nb_integration_pts);
1000 if (
dataAtPts->varGradRotAxis.size2() != nb_integration_pts)
1002 "varGradRotAxis columns (%d) != nb integration points (%d)",
1003 static_cast<int>(
dataAtPts->varGradRotAxis.size2()),
1004 nb_integration_pts);
1005 if (
dataAtPts->varPiola.size2() != nb_integration_pts)
1007 "varPiola columns (%d) != nb integration points (%d)",
1008 static_cast<int>(
dataAtPts->varPiola.size2()),
1009 nb_integration_pts);
1010 if (
dataAtPts->varDivPiola.size2() != nb_integration_pts)
1012 "varDivPiola columns (%d) != nb integration points (%d)",
1013 static_cast<int>(
dataAtPts->varDivPiola.size2()),
1014 nb_integration_pts);
1015#endif
1016
1019 auto t_div_P = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->divPAtPts);
1020 auto t_w_l2 = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->wL2AtPts);
1021 auto t_s_dot_w = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->wL2DotAtPts);
1022 auto t_s_dot_dot_w =
1023 getFTensor1FromMat<SPACE_DIM>(
dataAtPts->wL2DotDotAtPts);
1024 auto t_h = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->hAtPts);
1025 auto t_levi_kirchhoff =
1026 getFTensor1FromMat<SPACE_DIM>(
dataAtPts->leviKirchhoffAtPts);
1027 auto t_omega_grad_dot =
1028 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->rotAxisGradDotAtPts);
1029 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1030 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->stretchTensorAtPts);
1031
1032 auto t_var_w = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->varWL2);
1033 auto t_var_omega = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->varRotAxis);
1034 auto t_var_grad_omega =
1035 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->varGradRotAxis);
1036 auto t_var_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
dataAtPts->varPiola);
1037 auto t_var_div_P = getFTensor1FromMat<SPACE_DIM>(
dataAtPts->varDivPiola);
1038
1039 auto t_det = getFTensor0FromVec(
topoData->detJacobianAtPts);
1040 auto t_inv_jac =
1041 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
topoData->invJacobianAtPtr);
1042
1043 if (
dataAtPts->wL2DotDotAtPts.size1() == 0 &&
1044 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
1047 }
1048
1049 const auto piola_scale =
dataAtPts->piolaScale;
1050 const auto alpha_w =
alphaW / piola_scale;
1051 const auto alpha_rho =
alphaRho / piola_scale;
1052
1053 const int nb_base_functions = data.
getN().size2();
1055
1056 auto get_ftensor1 = [](
auto &
v) {
1058 &
v[0], &
v[1], &
v[2]);
1059 };
1060
1061 auto next = [&]() {
1062 ++t_w;
1063 ++t_div_P;
1064 ++t_w_l2;
1065 ++t_s_dot_w;
1066 ++t_s_dot_dot_w;
1067 ++t_h;
1068 ++t_levi_kirchhoff;
1069 ++t_omega_grad_dot;
1070 ++t_R;
1071 ++t_u;
1072
1073 ++t_var_w;
1074 ++t_var_omega;
1075 ++t_var_grad_omega;
1076 ++t_var_P;
1077 ++t_var_div_P;
1078
1079 ++t_det;
1080 ++t_inv_jac;
1081 };
1082
1089
1090 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1091
1092
1094 t_cof(
i,
j) = t_det * t_inv_jac(
j,
i);
1095
1096 auto t_nf = get_ftensor1(
nF);
1097 int bb = 0;
1098 for (; bb != nb_dofs /
SPACE_DIM; ++bb) {
1099
1101 t_div_base(
i) = -(1 / t_det) * (t_inv_jac(
j,
i) * t_base_diff(
j));
1102
1103
1104 t_nf(
i) += (t_w *
v) *
1105 (t_var_w(
k) * (-t_div_P(
k) + alpha_w * t_s_dot_w(
k) +
1106 alpha_rho * t_s_dot_dot_w(
k))) *
1107 t_cof(
i,
j) * t_base_diff(
j);
1108 t_nf(
i) += (t_w *
v) * (-(t_var_w(
k) * t_div_P(
k))) * t_div_base(
i);
1109
1110
1111 t_nf(
i) += (t_w *
v) * (t_var_omega(
k) * (-t_levi_kirchhoff(
k))) *
1112 t_cof(
i,
j) * t_base_diff(
j);
1113 #ifndef NDEBUG
1114
1116 SETERRQ(
1118 "OpSensitivity_dX with alpha_omega != 0 is not implemented yet");
1119 }
1120 #endif
1121
1122
1124 t_nf(
i) -= (t_w *
v) * (t_var_P(
i,
k) * (t_R(
i,
l) * t_u(
l,
k)) / 2) *
1125 t_cof(
i,
j) * t_base_diff(
j);
1126 t_nf(
i) -= (t_w *
v) * (t_var_P(
i,
l) * (t_R(
i,
k) * t_u(
l,
k)) / 2) *
1127 t_cof(
i,
j) * t_base_diff(
j);
1128 t_nf(
i) += (t_w *
v) * (t_var_P(
i,
j) *
t_kd(
i,
j)) * t_cof(
i,
j) *
1130 } else {
1133 t_nf(
i) += (t_w *
v) * (t_var_P(
k,
m) * (-t_residuum_P(
k,
m))) *
1134 t_cof(
i,
j) * t_base_diff(
j);
1135 }
1136
1137
1138 t_nf(
i) += (t_w *
v) * (t_var_div_P(
k) * (-t_w_l2(
k))) * t_cof(
i,
j) *
1140 t_nf(
i) += (t_w *
v) * (t_var_div_P(
k) * (-t_w_l2(
k))) * t_div_base(
i);
1141
1142 ++t_nf;
1143 ++t_base_diff;
1144 }
1145 for (; bb != nb_base_functions; ++bb)
1146 ++t_base_diff;
1147
1148 next();
1149 }
1150
1152 }
#define FTENSOR_INDEX(DIM, I)
#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 ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Index< 'm', 3 > m
static enum RotSelector gradApproximator
boost::shared_ptr< TopologicalData > topoData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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 degrees of freedom on entity.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getVolume() const
element volume (linear geometry)
VectorDouble nF
local right hand side vector
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts