668 namespace MappingFEFieldImplementation
678 template <
int dim,
int spacedim,
typename VectorType>
680 maybe_compute_q_points(
681 const typename ::QProjector<dim>::DataSetDescriptor data_set,
682 const typename ::MappingFEField<dim, spacedim, VectorType>::
686 const std::vector<unsigned int> &fe_to_real,
687 std::vector<Point<spacedim>> &quadrature_points)
693 for (
unsigned int point = 0; point < quadrature_points.size();
697 const double *shape = &data.shape(point + data_set, 0);
699 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
701 const unsigned int comp_k =
702 fe.system_to_component_index(k).first;
704 result[fe_to_real[comp_k]] +=
705 data.local_dof_values[k] * shape[k];
708 quadrature_points[point] = result;
720 template <
int dim,
int spacedim,
typename VectorType>
722 maybe_update_Jacobians(
723 const typename ::QProjector<dim>::DataSetDescriptor data_set,
724 const typename ::MappingFEField<dim, spacedim, VectorType>::
728 const std::vector<unsigned int> &fe_to_real)
735 const unsigned int n_q_points = data.contravariant.size();
739 for (
unsigned int point = 0; point < n_q_points; ++point)
742 &data.derivative(point + data_set, 0);
746 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
748 const unsigned int comp_k =
749 fe.system_to_component_index(k).first;
751 result[fe_to_real[comp_k]] +=
752 data.local_dof_values[k] * data_derv[k];
756 for (
unsigned int i = 0; i < spacedim; ++i)
758 data.contravariant[point][i] = result[i];
766 for (
unsigned int point = 0; point < data.contravariant.size();
768 data.covariant[point] =
769 (data.contravariant[point]).covariant_form();
775 data.volume_elements.size());
776 for (
unsigned int point = 0; point < data.contravariant.size();
778 data.volume_elements[point] =
779 data.contravariant[point].determinant();
789 template <
int dim,
int spacedim,
typename VectorType>
791 maybe_update_jacobian_grads(
792 const typename ::QProjector<dim>::DataSetDescriptor data_set,
793 const typename ::MappingFEField<dim, spacedim, VectorType>::
797 const std::vector<unsigned int> &fe_to_real,
798 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
803 const unsigned int n_q_points = jacobian_grads.size();
805 for (
unsigned int point = 0; point < n_q_points; ++point)
808 &data.second_derivative(point + data_set, 0);
812 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
814 const unsigned int comp_k =
815 fe.system_to_component_index(k).first;
817 for (
unsigned int j = 0; j < dim; ++j)
818 for (
unsigned int l = 0; l < dim; ++l)
819 result[fe_to_real[comp_k]][j][l] +=
820 (
second[k][j][l] * data.local_dof_values[k]);
825 for (
unsigned int i = 0; i < spacedim; ++i)
826 for (
unsigned int j = 0; j < dim; ++j)
827 for (
unsigned int l = 0; l < dim; ++l)
828 jacobian_grads[point][i][j][l] = result[i][j][l];
839 template <
int dim,
int spacedim,
typename VectorType>
841 maybe_update_jacobian_pushed_forward_grads(
842 const typename ::QProjector<dim>::DataSetDescriptor data_set,
843 const typename ::MappingFEField<dim, spacedim, VectorType>::
847 const std::vector<unsigned int> &fe_to_real,
848 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
853 const unsigned int n_q_points =
854 jacobian_pushed_forward_grads.size();
856 double tmp[spacedim][spacedim][spacedim];
857 for (
unsigned int point = 0; point < n_q_points; ++point)
860 &data.second_derivative(point + data_set, 0);
864 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
866 const unsigned int comp_k =
867 fe.system_to_component_index(k).first;
869 for (
unsigned int j = 0; j < dim; ++j)
870 for (
unsigned int l = 0; l < dim; ++l)
871 result[fe_to_real[comp_k]][j][l] +=
872 (
second[k][j][l] * data.local_dof_values[k]);
876 for (
unsigned int i = 0; i < spacedim; ++i)
877 for (
unsigned int j = 0; j < spacedim; ++j)
878 for (
unsigned int l = 0; l < dim; ++l)
881 result[i][0][l] * data.covariant[point][j][0];
882 for (
unsigned int jr = 1; jr < dim; ++jr)
885 result[i][jr][l] * data.covariant[point][j][jr];
890 for (
unsigned int i = 0; i < spacedim; ++i)
891 for (
unsigned int j = 0; j < spacedim; ++j)
892 for (
unsigned int l = 0; l < spacedim; ++l)
894 jacobian_pushed_forward_grads[point][i][j][l] =
895 tmp[i][j][0] * data.covariant[point][l][0];
896 for (
unsigned int lr = 1; lr < dim; ++lr)
898 jacobian_pushed_forward_grads[point][i][j][l] +=
899 tmp[i][j][lr] * data.covariant[point][l][lr];
912 template <
int dim,
int spacedim,
typename VectorType>
914 maybe_update_jacobian_2nd_derivatives(
915 const typename ::QProjector<dim>::DataSetDescriptor data_set,
916 const typename ::MappingFEField<dim, spacedim, VectorType>::
920 const std::vector<unsigned int> &fe_to_real,
921 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
926 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
928 for (
unsigned int point = 0; point < n_q_points; ++point)
931 &data.third_derivative(point + data_set, 0);
935 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
937 const unsigned int comp_k =
938 fe.system_to_component_index(k).first;
940 for (
unsigned int j = 0; j < dim; ++j)
941 for (
unsigned int l = 0; l < dim; ++l)
942 for (
unsigned int m = 0; m < dim; ++m)
943 result[fe_to_real[comp_k]][j][l][m] +=
944 (third[k][j][l][m] * data.local_dof_values[k]);
949 for (
unsigned int i = 0; i < spacedim; ++i)
950 for (
unsigned int j = 0; j < dim; ++j)
951 for (
unsigned int l = 0; l < dim; ++l)
952 for (
unsigned int m = 0; m < dim; ++m)
953 jacobian_2nd_derivatives[point][i][j][l][m] =
966 template <
int dim,
int spacedim,
typename VectorType>
968 maybe_update_jacobian_pushed_forward_2nd_derivatives(
969 const typename ::QProjector<dim>::DataSetDescriptor data_set,
970 const typename ::MappingFEField<dim, spacedim, VectorType>::
974 const std::vector<unsigned int> &fe_to_real,
975 std::vector<Tensor<4, spacedim>>
976 &jacobian_pushed_forward_2nd_derivatives)
981 const unsigned int n_q_points =
982 jacobian_pushed_forward_2nd_derivatives.size();
984 double tmp[spacedim][spacedim][spacedim][spacedim];
985 for (
unsigned int point = 0; point < n_q_points; ++point)
988 &data.third_derivative(point + data_set, 0);
992 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
994 const unsigned int comp_k =
995 fe.system_to_component_index(k).first;
997 for (
unsigned int j = 0; j < dim; ++j)
998 for (
unsigned int l = 0; l < dim; ++l)
999 for (
unsigned int m = 0; m < dim; ++m)
1000 result[fe_to_real[comp_k]][j][l][m] +=
1001 (third[k][j][l][m] * data.local_dof_values[k]);
1005 for (
unsigned int i = 0; i < spacedim; ++i)
1006 for (
unsigned int j = 0; j < spacedim; ++j)
1007 for (
unsigned int l = 0; l < dim; ++l)
1008 for (
unsigned int m = 0; m < dim; ++m)
1010 jacobian_pushed_forward_2nd_derivatives
1011 [point][i][j][l][m] =
1012 result[i][0][l][m] * data.covariant[point][j][0];
1013 for (
unsigned int jr = 1; jr < dim; ++jr)
1014 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1016 result[i][jr][l][m] *
1017 data.covariant[point][j][jr];
1021 for (
unsigned int i = 0; i < spacedim; ++i)
1022 for (
unsigned int j = 0; j < spacedim; ++j)
1023 for (
unsigned int l = 0; l < spacedim; ++l)
1024 for (
unsigned int m = 0; m < dim; ++m)
1027 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1029 data.covariant[point][l][0];
1030 for (
unsigned int lr = 1; lr < dim; ++lr)
1032 jacobian_pushed_forward_2nd_derivatives[point][i]
1035 data.covariant[point][l][lr];
1039 for (
unsigned int i = 0; i < spacedim; ++i)
1040 for (
unsigned int j = 0; j < spacedim; ++j)
1041 for (
unsigned int l = 0; l < spacedim; ++l)
1042 for (
unsigned int m = 0; m < spacedim; ++m)
1044 jacobian_pushed_forward_2nd_derivatives
1045 [point][i][j][l][m] =
1046 tmp[i][j][l][0] * data.covariant[point][m][0];
1047 for (
unsigned int mr = 1; mr < dim; ++mr)
1048 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1050 tmp[i][j][l][mr] * data.covariant[point][m][mr];
1062 template <
int dim,
int spacedim,
typename VectorType>
1064 maybe_update_jacobian_3rd_derivatives(
1065 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1066 const typename ::MappingFEField<dim, spacedim, VectorType>::
1070 const std::vector<unsigned int> &fe_to_real,
1071 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1073 const UpdateFlags update_flags = data.update_each;
1076 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1078 for (
unsigned int point = 0; point < n_q_points; ++point)
1081 &data.fourth_derivative(point + data_set, 0);
1085 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1087 const unsigned int comp_k =
1088 fe.system_to_component_index(k).first;
1089 if (fe_mask[comp_k])
1090 for (
unsigned int j = 0; j < dim; ++j)
1091 for (
unsigned int l = 0; l < dim; ++l)
1092 for (
unsigned int m = 0; m < dim; ++m)
1093 for (
unsigned int n = 0; n < dim; ++n)
1094 result[fe_to_real[comp_k]][j][l][m][n] +=
1095 (fourth[k][j][l][m][n] *
1096 data.local_dof_values[k]);
1102 for (
unsigned int i = 0; i < spacedim; ++i)
1103 for (
unsigned int j = 0; j < dim; ++j)
1104 for (
unsigned int l = 0; l < dim; ++l)
1105 for (
unsigned int m = 0; m < dim; ++m)
1106 for (
unsigned int n = 0; n < dim; ++n)
1107 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1108 result[i][j][l][m][n];
1120 template <
int dim,
int spacedim,
typename VectorType>
1122 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1123 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1124 const typename ::MappingFEField<dim, spacedim, VectorType>::
1128 const std::vector<unsigned int> &fe_to_real,
1129 std::vector<Tensor<5, spacedim>>
1130 &jacobian_pushed_forward_3rd_derivatives)
1132 const UpdateFlags update_flags = data.update_each;
1135 const unsigned int n_q_points =
1136 jacobian_pushed_forward_3rd_derivatives.size();
1138 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1139 for (
unsigned int point = 0; point < n_q_points; ++point)
1142 &data.fourth_derivative(point + data_set, 0);
1146 for (
unsigned int k = 0; k < data.n_shape_functions; ++k)
1148 const unsigned int comp_k =
1149 fe.system_to_component_index(k).first;
1150 if (fe_mask[comp_k])
1151 for (
unsigned int j = 0; j < dim; ++j)
1152 for (
unsigned int l = 0; l < dim; ++l)
1153 for (
unsigned int m = 0; m < dim; ++m)
1154 for (
unsigned int n = 0; n < dim; ++n)
1155 result[fe_to_real[comp_k]][j][l][m][n] +=
1156 (fourth[k][j][l][m][n] *
1157 data.local_dof_values[k]);
1161 for (
unsigned int i = 0; i < spacedim; ++i)
1162 for (
unsigned int j = 0; j < spacedim; ++j)
1163 for (
unsigned int l = 0; l < dim; ++l)
1164 for (
unsigned int m = 0; m < dim; ++m)
1165 for (
unsigned int n = 0; n < dim; ++n)
1167 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1168 data.covariant[point][j][0];
1169 for (
unsigned int jr = 1; jr < dim; ++jr)
1170 tmp[i][j][l][m][n] +=
1171 result[i][jr][l][m][n] *
1172 data.covariant[point][j][jr];
1176 for (
unsigned int i = 0; i < spacedim; ++i)
1177 for (
unsigned int j = 0; j < spacedim; ++j)
1178 for (
unsigned int l = 0; l < spacedim; ++l)
1179 for (
unsigned int m = 0; m < dim; ++m)
1180 for (
unsigned int n = 0; n < dim; ++n)
1182 jacobian_pushed_forward_3rd_derivatives
1183 [point][i][j][l][m][n] =
1184 tmp[i][j][0][m][n] *
1185 data.covariant[point][l][0];
1186 for (
unsigned int lr = 1; lr < dim; ++lr)
1187 jacobian_pushed_forward_3rd_derivatives[point][i]
1190 tmp[i][j][lr][m][n] *
1191 data.covariant[point][l][lr];
1195 for (
unsigned int i = 0; i < spacedim; ++i)
1196 for (
unsigned int j = 0; j < spacedim; ++j)
1197 for (
unsigned int l = 0; l < spacedim; ++l)
1198 for (
unsigned int m = 0; m < spacedim; ++m)
1199 for (
unsigned int n = 0; n < dim; ++n)
1201 tmp[i][j][l][m][n] =
1202 jacobian_pushed_forward_3rd_derivatives[point][i]
1205 data.covariant[point][m][0];
1206 for (
unsigned int mr = 1; mr < dim; ++mr)
1207 tmp[i][j][l][m][n] +=
1208 jacobian_pushed_forward_3rd_derivatives[point]
1211 data.covariant[point][m][mr];
1215 for (
unsigned int i = 0; i < spacedim; ++i)
1216 for (
unsigned int j = 0; j < spacedim; ++j)
1217 for (
unsigned int l = 0; l < spacedim; ++l)
1218 for (
unsigned int m = 0; m < spacedim; ++m)
1219 for (
unsigned int n = 0; n < spacedim; ++n)
1221 jacobian_pushed_forward_3rd_derivatives
1222 [point][i][j][l][m][n] =
1223 tmp[i][j][l][m][0] *
1224 data.covariant[point][n][0];
1225 for (
unsigned int nr = 1; nr < dim; ++nr)
1226 jacobian_pushed_forward_3rd_derivatives[point][i]
1229 tmp[i][j][l][m][nr] *
1230 data.covariant[point][n][nr];
1246 template <
int dim,
int spacedim,
typename VectorType>
1248 maybe_compute_face_data(
1249 const ::Mapping<dim, spacedim> &mapping,
1250 const typename ::Triangulation<dim, spacedim>::cell_iterator
1252 const unsigned int face_no,
1253 const unsigned int subface_no,
1255 const typename ::MappingFEField<dim, spacedim, VectorType>::
1260 const UpdateFlags update_flags = data.update_each;
1264 const unsigned int n_q_points = output_data.boundary_forms.size();
1273 for (
unsigned int d = 0; d != dim - 1; ++d)
1275 Assert(face_no + cell->n_faces() * d <
1276 data.unit_tangentials.size(),
1279 data.aux[d].size() <=
1280 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1285 data.unit_tangentials[face_no + cell->n_faces() * d]),
1293 if (dim == spacedim)
1295 for (
unsigned int i = 0; i < n_q_points; ++i)
1303 output_data.boundary_forms[i][0] =
1304 (face_no == 0 ? -1 : +1);
1307 output_data.boundary_forms[i] =
1308 cross_product_2d(data.aux[0][i]);
1311 output_data.boundary_forms[i] =
1312 cross_product_3d(data.aux[0][i], data.aux[1][i]);
1328 for (
unsigned int point = 0; point < n_q_points; ++point)
1333 output_data.boundary_forms[point] =
1334 data.contravariant[point].transpose()[0];
1335 output_data.boundary_forms[point] /=
1336 (face_no == 0 ? -1. : +1.) *
1337 output_data.boundary_forms[point].norm();
1346 cross_product_3d(DX_t[0], DX_t[1]);
1347 cell_normal /= cell_normal.
norm();
1351 output_data.boundary_forms[point] =
1352 cross_product_3d(data.aux[0][point], cell_normal);
1358 for (
unsigned int i = 0; i < output_data.boundary_forms.size();
1363 output_data.JxW_values[i] =
1364 output_data.boundary_forms[i].norm() *
1365 data.quadrature_weights[i + data_set];
1370 const double area_ratio =
1372 cell->subface_case(face_no), subface_no);
1373 output_data.JxW_values[i] *= area_ratio;
1378 output_data.normal_vectors[i] =
1380 output_data.boundary_forms[i].norm());
1391 template <
int dim,
int spacedim,
typename VectorType>
1393 do_fill_fe_face_values(
1394 const ::Mapping<dim, spacedim> &mapping,
1395 const typename ::Triangulation<dim, spacedim>::cell_iterator
1397 const unsigned int face_no,
1398 const unsigned int subface_no,
1399 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1400 const typename ::MappingFEField<dim, spacedim, VectorType>::
1404 const std::vector<unsigned int> &fe_to_real,
1408 maybe_compute_q_points<dim, spacedim, VectorType>(
1414 output_data.quadrature_points);
1416 maybe_update_Jacobians<dim, spacedim, VectorType>(
1417 data_set, data, fe, fe_mask, fe_to_real);
1419 const UpdateFlags update_flags = data.update_each;
1420 const unsigned int n_q_points = data.contravariant.size();
1423 for (
unsigned int point = 0; point < n_q_points; ++point)
1424 output_data.jacobians[point] = data.contravariant[point];
1427 for (
unsigned int point = 0; point < n_q_points; ++point)
1428 output_data.inverse_jacobians[point] =
1429 data.covariant[point].transpose();
1431 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1432 data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1434 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1440 output_data.jacobian_pushed_forward_grads);
1442 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1448 output_data.jacobian_2nd_derivatives);
1450 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1458 output_data.jacobian_pushed_forward_2nd_derivatives);
1460 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1466 output_data.jacobian_3rd_derivatives);
1468 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1476 output_data.jacobian_pushed_forward_3rd_derivatives);
1478 maybe_compute_face_data<dim, spacedim, VectorType>(
1479 mapping, cell, face_no, subface_no, data_set, data, output_data);