671 namespace MappingFEFieldImplementation
681 template <
int dim,
int spacedim,
typename VectorType>
683 maybe_compute_q_points(
684 const typename ::QProjector<dim>::DataSetDescriptor data_set,
685 const typename ::MappingFEField<dim, spacedim, VectorType>::
689 const std::vector<unsigned int> & fe_to_real,
690 std::vector<Point<spacedim>> & quadrature_points)
696 for (
unsigned int point = 0; point < quadrature_points.size();
700 const double * shape = &data.shape(point + data_set, 0);
702 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
704 const unsigned int comp_k =
705 fe.system_to_component_index(k).first;
707 result[fe_to_real[comp_k]] +=
708 data.local_dof_values[k] * shape[k];
711 quadrature_points[point] = result;
723 template <
int dim,
int spacedim,
typename VectorType>
725 maybe_update_Jacobians(
726 const typename ::QProjector<dim>::DataSetDescriptor data_set,
727 const typename ::MappingFEField<dim, spacedim, VectorType>::
731 const std::vector<unsigned int> & fe_to_real)
738 const unsigned int n_q_points = data.contravariant.size();
742 for (
unsigned int point = 0; point < n_q_points; ++point)
745 &data.derivative(point + data_set, 0);
749 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
751 const unsigned int comp_k =
752 fe.system_to_component_index(k).first;
754 result[fe_to_real[comp_k]] +=
755 data.local_dof_values[k] * data_derv[k];
759 for (
unsigned int i = 0; i < spacedim; ++i)
761 data.contravariant[point][i] = result[i];
769 for (
unsigned int point = 0; point < data.contravariant.size();
771 data.covariant[point] =
772 (data.contravariant[point]).covariant_form();
778 data.volume_elements.size());
779 for (
unsigned int point = 0; point < data.contravariant.size();
781 data.volume_elements[point] =
782 data.contravariant[point].determinant();
792 template <
int dim,
int spacedim,
typename VectorType>
794 maybe_update_jacobian_grads(
795 const typename ::QProjector<dim>::DataSetDescriptor data_set,
796 const typename ::MappingFEField<dim, spacedim, VectorType>::
800 const std::vector<unsigned int> & fe_to_real,
801 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
806 const unsigned int n_q_points = jacobian_grads.size();
808 for (
unsigned int point = 0; point < n_q_points; ++point)
811 &data.second_derivative(point + data_set, 0);
815 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
817 const unsigned int comp_k =
818 fe.system_to_component_index(k).first;
820 for (
unsigned int j = 0; j < dim; ++j)
821 for (
unsigned int l = 0; l < dim; ++l)
822 result[fe_to_real[comp_k]][j][l] +=
823 (
second[k][j][l] * data.local_dof_values[k]);
828 for (
unsigned int i = 0; i < spacedim; ++i)
829 for (
unsigned int j = 0; j < dim; ++j)
830 for (
unsigned int l = 0; l < dim; ++l)
831 jacobian_grads[point][i][j][l] = result[i][j][l];
842 template <
int dim,
int spacedim,
typename VectorType>
844 maybe_update_jacobian_pushed_forward_grads(
845 const typename ::QProjector<dim>::DataSetDescriptor data_set,
846 const typename ::MappingFEField<dim, spacedim, VectorType>::
850 const std::vector<unsigned int> & fe_to_real,
851 std::vector<Tensor<3, spacedim>> & jacobian_pushed_forward_grads)
856 const unsigned int n_q_points =
857 jacobian_pushed_forward_grads.size();
859 double tmp[spacedim][spacedim][spacedim];
860 for (
unsigned int point = 0; point < n_q_points; ++point)
863 &data.second_derivative(point + data_set, 0);
867 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
869 const unsigned int comp_k =
870 fe.system_to_component_index(k).first;
872 for (
unsigned int j = 0; j < dim; ++j)
873 for (
unsigned int l = 0; l < dim; ++l)
874 result[fe_to_real[comp_k]][j][l] +=
875 (
second[k][j][l] * data.local_dof_values[k]);
879 for (
unsigned int i = 0; i < spacedim; ++i)
880 for (
unsigned int j = 0; j < spacedim; ++j)
881 for (
unsigned int l = 0; l < dim; ++l)
884 result[i][0][l] * data.covariant[point][j][0];
885 for (
unsigned int jr = 1; jr < dim; ++jr)
888 result[i][jr][l] * data.covariant[point][j][jr];
893 for (
unsigned int i = 0; i < spacedim; ++i)
894 for (
unsigned int j = 0; j < spacedim; ++j)
895 for (
unsigned int l = 0; l < spacedim; ++l)
897 jacobian_pushed_forward_grads[point][i][j][l] =
898 tmp[i][j][0] * data.covariant[point][l][0];
899 for (
unsigned int lr = 1; lr < dim; ++lr)
901 jacobian_pushed_forward_grads[point][i][j][l] +=
902 tmp[i][j][lr] * data.covariant[point][l][lr];
915 template <
int dim,
int spacedim,
typename VectorType>
917 maybe_update_jacobian_2nd_derivatives(
918 const typename ::QProjector<dim>::DataSetDescriptor data_set,
919 const typename ::MappingFEField<dim, spacedim, VectorType>::
923 const std::vector<unsigned int> & fe_to_real,
924 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
929 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
931 for (
unsigned int point = 0; point < n_q_points; ++point)
934 &data.third_derivative(point + data_set, 0);
938 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
940 const unsigned int comp_k =
941 fe.system_to_component_index(k).first;
943 for (
unsigned int j = 0; j < dim; ++j)
944 for (
unsigned int l = 0; l < dim; ++l)
945 for (
unsigned int m = 0; m < dim; ++m)
946 result[fe_to_real[comp_k]][j][l][m] +=
947 (third[k][j][l][m] * data.local_dof_values[k]);
952 for (
unsigned int i = 0; i < spacedim; ++i)
953 for (
unsigned int j = 0; j < dim; ++j)
954 for (
unsigned int l = 0; l < dim; ++l)
955 for (
unsigned int m = 0; m < dim; ++m)
956 jacobian_2nd_derivatives[point][i][j][l][m] =
969 template <
int dim,
int spacedim,
typename VectorType>
971 maybe_update_jacobian_pushed_forward_2nd_derivatives(
972 const typename ::QProjector<dim>::DataSetDescriptor data_set,
973 const typename ::MappingFEField<dim, spacedim, VectorType>::
977 const std::vector<unsigned int> & fe_to_real,
978 std::vector<Tensor<4, spacedim>>
979 &jacobian_pushed_forward_2nd_derivatives)
984 const unsigned int n_q_points =
985 jacobian_pushed_forward_2nd_derivatives.size();
987 double tmp[spacedim][spacedim][spacedim][spacedim];
988 for (
unsigned int point = 0; point < n_q_points; ++point)
991 &data.third_derivative(point + data_set, 0);
995 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
997 const unsigned int comp_k =
998 fe.system_to_component_index(k).first;
1000 for (
unsigned int j = 0; j < dim; ++j)
1001 for (
unsigned int l = 0; l < dim; ++l)
1002 for (
unsigned int m = 0; m < dim; ++m)
1003 result[fe_to_real[comp_k]][j][l][m] +=
1004 (third[k][j][l][m] * data.local_dof_values[k]);
1008 for (
unsigned int i = 0; i < spacedim; ++i)
1009 for (
unsigned int j = 0; j < spacedim; ++j)
1010 for (
unsigned int l = 0; l < dim; ++l)
1011 for (
unsigned int m = 0; m < dim; ++m)
1013 jacobian_pushed_forward_2nd_derivatives
1014 [point][i][j][l][m] =
1015 result[i][0][l][m] * data.covariant[point][j][0];
1016 for (
unsigned int jr = 1; jr < dim; ++jr)
1017 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1019 result[i][jr][l][m] *
1020 data.covariant[point][j][jr];
1024 for (
unsigned int i = 0; i < spacedim; ++i)
1025 for (
unsigned int j = 0; j < spacedim; ++j)
1026 for (
unsigned int l = 0; l < spacedim; ++l)
1027 for (
unsigned int m = 0; m < dim; ++m)
1030 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1032 data.covariant[point][l][0];
1033 for (
unsigned int lr = 1; lr < dim; ++lr)
1035 jacobian_pushed_forward_2nd_derivatives[point][i]
1038 data.covariant[point][l][lr];
1042 for (
unsigned int i = 0; i < spacedim; ++i)
1043 for (
unsigned int j = 0; j < spacedim; ++j)
1044 for (
unsigned int l = 0; l < spacedim; ++l)
1045 for (
unsigned int m = 0; m < spacedim; ++m)
1047 jacobian_pushed_forward_2nd_derivatives
1048 [point][i][j][l][m] =
1049 tmp[i][j][l][0] * data.covariant[point][m][0];
1050 for (
unsigned int mr = 1; mr < dim; ++mr)
1051 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1053 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1065 template <
int dim,
int spacedim,
typename VectorType>
1067 maybe_update_jacobian_3rd_derivatives(
1068 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1069 const typename ::MappingFEField<dim, spacedim, VectorType>::
1070 InternalData & data,
1073 const std::vector<unsigned int> & fe_to_real,
1074 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1076 const UpdateFlags update_flags = data.update_each;
1079 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1081 for (
unsigned int point = 0; point < n_q_points; ++point)
1084 &data.fourth_derivative(point + data_set, 0);
1088 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1090 const unsigned int comp_k =
1091 fe.system_to_component_index(k).first;
1092 if (fe_mask[comp_k])
1093 for (
unsigned int j = 0; j < dim; ++j)
1094 for (
unsigned int l = 0; l < dim; ++l)
1095 for (
unsigned int m = 0; m < dim; ++m)
1096 for (
unsigned int n = 0; n < dim; ++n)
1097 result[fe_to_real[comp_k]][j][l][m][n] +=
1098 (fourth[k][j][l][m][n] *
1099 data.local_dof_values[k]);
1105 for (
unsigned int i = 0; i < spacedim; ++i)
1106 for (
unsigned int j = 0; j < dim; ++j)
1107 for (
unsigned int l = 0; l < dim; ++l)
1108 for (
unsigned int m = 0; m < dim; ++m)
1109 for (
unsigned int n = 0; n < dim; ++n)
1110 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1111 result[i][j][l][m][n];
1123 template <
int dim,
int spacedim,
typename VectorType>
1125 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1126 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1127 const typename ::MappingFEField<dim, spacedim, VectorType>::
1128 InternalData & data,
1131 const std::vector<unsigned int> & fe_to_real,
1132 std::vector<Tensor<5, spacedim>>
1133 &jacobian_pushed_forward_3rd_derivatives)
1135 const UpdateFlags update_flags = data.update_each;
1138 const unsigned int n_q_points =
1139 jacobian_pushed_forward_3rd_derivatives.size();
1141 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1142 for (
unsigned int point = 0; point < n_q_points; ++point)
1145 &data.fourth_derivative(point + data_set, 0);
1149 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1151 const unsigned int comp_k =
1152 fe.system_to_component_index(k).first;
1153 if (fe_mask[comp_k])
1154 for (
unsigned int j = 0; j < dim; ++j)
1155 for (
unsigned int l = 0; l < dim; ++l)
1156 for (
unsigned int m = 0; m < dim; ++m)
1157 for (
unsigned int n = 0; n < dim; ++n)
1158 result[fe_to_real[comp_k]][j][l][m][n] +=
1159 (fourth[k][j][l][m][n] *
1160 data.local_dof_values[k]);
1164 for (
unsigned int i = 0; i < spacedim; ++i)
1165 for (
unsigned int j = 0; j < spacedim; ++j)
1166 for (
unsigned int l = 0; l < dim; ++l)
1167 for (
unsigned int m = 0; m < dim; ++m)
1168 for (
unsigned int n = 0; n < dim; ++n)
1170 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1171 data.covariant[point][j][0];
1172 for (
unsigned int jr = 1; jr < dim; ++jr)
1173 tmp[i][j][l][m][n] +=
1174 result[i][jr][l][m][n] *
1175 data.covariant[point][j][jr];
1179 for (
unsigned int i = 0; i < spacedim; ++i)
1180 for (
unsigned int j = 0; j < spacedim; ++j)
1181 for (
unsigned int l = 0; l < spacedim; ++l)
1182 for (
unsigned int m = 0; m < dim; ++m)
1183 for (
unsigned int n = 0; n < dim; ++n)
1185 jacobian_pushed_forward_3rd_derivatives
1186 [point][i][j][l][m][n] =
1187 tmp[i][j][0][m][n] *
1188 data.covariant[point][l][0];
1189 for (
unsigned int lr = 1; lr < dim; ++lr)
1190 jacobian_pushed_forward_3rd_derivatives[point][i]
1193 tmp[i][j][lr][m][n] *
1194 data.covariant[point][l][lr];
1198 for (
unsigned int i = 0; i < spacedim; ++i)
1199 for (
unsigned int j = 0; j < spacedim; ++j)
1200 for (
unsigned int l = 0; l < spacedim; ++l)
1201 for (
unsigned int m = 0; m < spacedim; ++m)
1202 for (
unsigned int n = 0; n < dim; ++n)
1204 tmp[i][j][l][m][n] =
1205 jacobian_pushed_forward_3rd_derivatives[point][i]
1208 data.covariant[point][m][0];
1209 for (
unsigned int mr = 1; mr < dim; ++mr)
1210 tmp[i][j][l][m][n] +=
1211 jacobian_pushed_forward_3rd_derivatives[point]
1214 data.covariant[point][m][mr];
1218 for (
unsigned int i = 0; i < spacedim; ++i)
1219 for (
unsigned int j = 0; j < spacedim; ++j)
1220 for (
unsigned int l = 0; l < spacedim; ++l)
1221 for (
unsigned int m = 0; m < spacedim; ++m)
1222 for (
unsigned int n = 0; n < spacedim; ++n)
1224 jacobian_pushed_forward_3rd_derivatives
1225 [point][i][j][l][m][n] =
1226 tmp[i][j][l][m][0] *
1227 data.covariant[point][n][0];
1228 for (
unsigned int nr = 1; nr < dim; ++nr)
1229 jacobian_pushed_forward_3rd_derivatives[point][i]
1232 tmp[i][j][l][m][nr] *
1233 data.covariant[point][n][nr];
1249 template <
int dim,
int spacedim,
typename VectorType>
1251 maybe_compute_face_data(
1252 const ::Mapping<dim, spacedim> &mapping,
1253 const typename ::Triangulation<dim, spacedim>::cell_iterator
1255 const unsigned int face_no,
1256 const unsigned int subface_no,
1258 const typename ::MappingFEField<dim, spacedim, VectorType>::
1263 const UpdateFlags update_flags = data.update_each;
1267 const unsigned int n_q_points = output_data.boundary_forms.size();
1276 for (
unsigned int d = 0; d != dim - 1; ++d)
1278 Assert(face_no + cell->n_faces() * d <
1279 data.unit_tangentials.size(),
1282 data.aux[d].size() <=
1283 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1288 data.unit_tangentials[face_no + cell->n_faces() * d]),
1296 if (dim == spacedim)
1298 for (
unsigned int i = 0; i < n_q_points; ++i)
1306 output_data.boundary_forms[i][0] =
1307 (face_no == 0 ? -1 : +1);
1310 output_data.boundary_forms[i] =
1311 cross_product_2d(data.aux[0][i]);
1314 output_data.boundary_forms[i] =
1315 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1331 for (
unsigned int point = 0; point < n_q_points; ++point)
1336 output_data.boundary_forms[point] =
1337 data.contravariant[point].transpose()[0];
1338 output_data.boundary_forms[point] /=
1339 (face_no == 0 ? -1. : +1.) *
1340 output_data.boundary_forms[point].norm();
1349 cross_product_3d(DX_t[0], DX_t[1]);
1350 cell_normal /= cell_normal.
norm();
1354 output_data.boundary_forms[point] =
1355 cross_product_3d(data.aux[0][point], cell_normal);
1361 for (
unsigned int i = 0; i < output_data.boundary_forms.size();
1366 output_data.JxW_values[i] =
1367 output_data.boundary_forms[i].norm() *
1368 data.quadrature_weights[i + data_set];
1373 const double area_ratio =
1375 cell->subface_case(face_no), subface_no);
1376 output_data.JxW_values[i] *= area_ratio;
1381 output_data.normal_vectors[i] =
1383 output_data.boundary_forms[i].norm());
1394 template <
int dim,
int spacedim,
typename VectorType>
1396 do_fill_fe_face_values(
1397 const ::Mapping<dim, spacedim> &mapping,
1398 const typename ::Triangulation<dim, spacedim>::cell_iterator
1400 const unsigned int face_no,
1401 const unsigned int subface_no,
1402 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1403 const typename ::MappingFEField<dim, spacedim, VectorType>::
1404 InternalData & data,
1407 const std::vector<unsigned int> & fe_to_real,
1411 maybe_compute_q_points<dim, spacedim, VectorType>(
1417 output_data.quadrature_points);
1419 maybe_update_Jacobians<dim, spacedim, VectorType>(
1420 data_set, data, fe, fe_mask, fe_to_real);
1422 const UpdateFlags update_flags = data.update_each;
1423 const unsigned int n_q_points = data.contravariant.size();
1426 for (
unsigned int point = 0; point < n_q_points; ++point)
1427 output_data.jacobians[point] = data.contravariant[point];
1430 for (
unsigned int point = 0; point < n_q_points; ++point)
1431 output_data.inverse_jacobians[point] =
1432 data.covariant[point].transpose();
1434 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1435 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1437 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1443 output_data.jacobian_pushed_forward_grads);
1445 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1451 output_data.jacobian_2nd_derivatives);
1453 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1461 output_data.jacobian_pushed_forward_2nd_derivatives);
1463 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1469 output_data.jacobian_3rd_derivatives);
1471 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1479 output_data.jacobian_pushed_forward_3rd_derivatives);
1481 maybe_compute_face_data<dim, spacedim, VectorType>(
1482 mapping, cell, face_no, subface_no, data_set, data, output_data);