684 namespace MappingFEFieldImplementation
694 template <
int dim,
int spacedim,
typename VectorType>
696 maybe_compute_q_points(
697 const typename ::QProjector<dim>::DataSetDescriptor data_set,
698 const typename ::MappingFEField<dim, spacedim, VectorType>::
702 const std::vector<unsigned int> &fe_to_real,
703 std::vector<Point<spacedim>> &quadrature_points)
709 for (
unsigned int point = 0; point < quadrature_points.size();
713 const double *shape = &
data.shape(point + data_set, 0);
715 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
717 const unsigned int comp_k =
718 fe.system_to_component_index(k).first;
720 result[fe_to_real[comp_k]] +=
721 data.local_dof_values[k] * shape[k];
724 quadrature_points[point] = result;
736 template <
int dim,
int spacedim,
typename VectorType>
738 maybe_update_Jacobians(
739 const typename ::QProjector<dim>::DataSetDescriptor data_set,
740 const typename ::MappingFEField<dim, spacedim, VectorType>::
744 const std::vector<unsigned int> &fe_to_real)
751 const unsigned int n_q_points =
data.contravariant.size();
755 for (
unsigned int point = 0; point < n_q_points; ++point)
758 &
data.derivative(point + data_set, 0);
762 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
764 const unsigned int comp_k =
765 fe.system_to_component_index(k).first;
767 result[fe_to_real[comp_k]] +=
768 data.local_dof_values[k] * data_derv[k];
772 for (
unsigned int i = 0; i < spacedim; ++i)
774 data.contravariant[point][i] = result[i];
782 for (
unsigned int point = 0; point <
data.contravariant.size();
784 data.covariant[point] =
785 (
data.contravariant[point]).covariant_form();
791 data.volume_elements.size());
792 for (
unsigned int point = 0; point <
data.contravariant.size();
794 data.volume_elements[point] =
795 data.contravariant[point].determinant();
805 template <
int dim,
int spacedim,
typename VectorType>
807 maybe_update_jacobian_grads(
808 const typename ::QProjector<dim>::DataSetDescriptor data_set,
809 const typename ::MappingFEField<dim, spacedim, VectorType>::
813 const std::vector<unsigned int> &fe_to_real,
814 std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
819 const unsigned int n_q_points = jacobian_grads.size();
821 for (
unsigned int point = 0; point < n_q_points; ++point)
824 &
data.second_derivative(point + data_set, 0);
828 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
830 const unsigned int comp_k =
831 fe.system_to_component_index(k).first;
833 for (
unsigned int j = 0; j < dim; ++j)
834 for (
unsigned int l = 0; l < dim; ++l)
835 result[fe_to_real[comp_k]][j][l] +=
841 for (
unsigned int i = 0; i < spacedim; ++i)
842 for (
unsigned int j = 0; j < dim; ++j)
843 for (
unsigned int l = 0; l < dim; ++l)
844 jacobian_grads[point][i][j][l] = result[i][j][l];
855 template <
int dim,
int spacedim,
typename VectorType>
857 maybe_update_jacobian_pushed_forward_grads(
858 const typename ::QProjector<dim>::DataSetDescriptor data_set,
859 const typename ::MappingFEField<dim, spacedim, VectorType>::
863 const std::vector<unsigned int> &fe_to_real,
864 std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
869 const unsigned int n_q_points =
870 jacobian_pushed_forward_grads.size();
872 double tmp[spacedim][spacedim][spacedim];
873 for (
unsigned int point = 0; point < n_q_points; ++point)
876 &
data.second_derivative(point + data_set, 0);
880 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
882 const unsigned int comp_k =
883 fe.system_to_component_index(k).first;
885 for (
unsigned int j = 0; j < dim; ++j)
886 for (
unsigned int l = 0; l < dim; ++l)
887 result[fe_to_real[comp_k]][j][l] +=
892 for (
unsigned int i = 0; i < spacedim; ++i)
893 for (
unsigned int j = 0; j < spacedim; ++j)
894 for (
unsigned int l = 0; l < dim; ++l)
897 result[i][0][l] *
data.covariant[point][j][0];
898 for (
unsigned int jr = 1; jr < dim; ++jr)
901 result[i][jr][l] *
data.covariant[point][j][jr];
906 for (
unsigned int i = 0; i < spacedim; ++i)
907 for (
unsigned int j = 0; j < spacedim; ++j)
908 for (
unsigned int l = 0; l < spacedim; ++l)
910 jacobian_pushed_forward_grads[point][i][j][l] =
911 tmp[i][j][0] *
data.covariant[point][l][0];
912 for (
unsigned int lr = 1; lr < dim; ++lr)
914 jacobian_pushed_forward_grads[point][i][j][l] +=
915 tmp[i][j][lr] *
data.covariant[point][l][lr];
928 template <
int dim,
int spacedim,
typename VectorType>
930 maybe_update_jacobian_2nd_derivatives(
931 const typename ::QProjector<dim>::DataSetDescriptor data_set,
932 const typename ::MappingFEField<dim, spacedim, VectorType>::
936 const std::vector<unsigned int> &fe_to_real,
937 std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
942 const unsigned int n_q_points = jacobian_2nd_derivatives.size();
944 for (
unsigned int point = 0; point < n_q_points; ++point)
947 &
data.third_derivative(point + data_set, 0);
951 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
953 const unsigned int comp_k =
954 fe.system_to_component_index(k).first;
956 for (
unsigned int j = 0; j < dim; ++j)
957 for (
unsigned int l = 0; l < dim; ++l)
958 for (
unsigned int m = 0; m < dim; ++m)
959 result[fe_to_real[comp_k]][j][l][m] +=
960 (third[k][j][l][m] *
data.local_dof_values[k]);
965 for (
unsigned int i = 0; i < spacedim; ++i)
966 for (
unsigned int j = 0; j < dim; ++j)
967 for (
unsigned int l = 0; l < dim; ++l)
968 for (
unsigned int m = 0; m < dim; ++m)
969 jacobian_2nd_derivatives[point][i][j][l][m] =
982 template <
int dim,
int spacedim,
typename VectorType>
984 maybe_update_jacobian_pushed_forward_2nd_derivatives(
985 const typename ::QProjector<dim>::DataSetDescriptor data_set,
986 const typename ::MappingFEField<dim, spacedim, VectorType>::
990 const std::vector<unsigned int> &fe_to_real,
991 std::vector<Tensor<4, spacedim>>
992 &jacobian_pushed_forward_2nd_derivatives)
997 const unsigned int n_q_points =
998 jacobian_pushed_forward_2nd_derivatives.size();
1000 double tmp[spacedim][spacedim][spacedim][spacedim];
1001 for (
unsigned int point = 0; point < n_q_points; ++point)
1004 &
data.third_derivative(point + data_set, 0);
1008 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
1010 const unsigned int comp_k =
1011 fe.system_to_component_index(k).first;
1012 if (fe_mask[comp_k])
1013 for (
unsigned int j = 0; j < dim; ++j)
1014 for (
unsigned int l = 0; l < dim; ++l)
1015 for (
unsigned int m = 0; m < dim; ++m)
1016 result[fe_to_real[comp_k]][j][l][m] +=
1017 (third[k][j][l][m] *
data.local_dof_values[k]);
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 < dim; ++l)
1024 for (
unsigned int m = 0; m < dim; ++m)
1026 jacobian_pushed_forward_2nd_derivatives
1027 [point][i][j][l][m] =
1028 result[i][0][l][m] *
data.covariant[point][j][0];
1029 for (
unsigned int jr = 1; jr < dim; ++jr)
1030 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1032 result[i][jr][l][m] *
1033 data.covariant[point][j][jr];
1037 for (
unsigned int i = 0; i < spacedim; ++i)
1038 for (
unsigned int j = 0; j < spacedim; ++j)
1039 for (
unsigned int l = 0; l < spacedim; ++l)
1040 for (
unsigned int m = 0; m < dim; ++m)
1043 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1045 data.covariant[point][l][0];
1046 for (
unsigned int lr = 1; lr < dim; ++lr)
1048 jacobian_pushed_forward_2nd_derivatives[point][i]
1051 data.covariant[point][l][lr];
1055 for (
unsigned int i = 0; i < spacedim; ++i)
1056 for (
unsigned int j = 0; j < spacedim; ++j)
1057 for (
unsigned int l = 0; l < spacedim; ++l)
1058 for (
unsigned int m = 0; m < spacedim; ++m)
1060 jacobian_pushed_forward_2nd_derivatives
1061 [point][i][j][l][m] =
1062 tmp[i][j][l][0] *
data.covariant[point][m][0];
1063 for (
unsigned int mr = 1; mr < dim; ++mr)
1064 jacobian_pushed_forward_2nd_derivatives[point][i][j]
1066 tmp[i][j][l][mr] *
data.covariant[point][m][mr];
1078 template <
int dim,
int spacedim,
typename VectorType>
1080 maybe_update_jacobian_3rd_derivatives(
1081 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1082 const typename ::MappingFEField<dim, spacedim, VectorType>::
1086 const std::vector<unsigned int> &fe_to_real,
1087 std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
1092 const unsigned int n_q_points = jacobian_3rd_derivatives.size();
1094 for (
unsigned int point = 0; point < n_q_points; ++point)
1097 &
data.fourth_derivative(point + data_set, 0);
1101 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
1103 const unsigned int comp_k =
1104 fe.system_to_component_index(k).first;
1105 if (fe_mask[comp_k])
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 result[fe_to_real[comp_k]][j][l][m][n] +=
1111 (fourth[k][j][l][m][n] *
1112 data.local_dof_values[k]);
1118 for (
unsigned int i = 0; i < spacedim; ++i)
1119 for (
unsigned int j = 0; j < dim; ++j)
1120 for (
unsigned int l = 0; l < dim; ++l)
1121 for (
unsigned int m = 0; m < dim; ++m)
1122 for (
unsigned int n = 0; n < dim; ++n)
1123 jacobian_3rd_derivatives[point][i][j][l][m][n] =
1124 result[i][j][l][m][n];
1136 template <
int dim,
int spacedim,
typename VectorType>
1138 maybe_update_jacobian_pushed_forward_3rd_derivatives(
1139 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1140 const typename ::MappingFEField<dim, spacedim, VectorType>::
1144 const std::vector<unsigned int> &fe_to_real,
1145 std::vector<Tensor<5, spacedim>>
1146 &jacobian_pushed_forward_3rd_derivatives)
1151 const unsigned int n_q_points =
1152 jacobian_pushed_forward_3rd_derivatives.size();
1154 double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
1155 for (
unsigned int point = 0; point < n_q_points; ++point)
1158 &
data.fourth_derivative(point + data_set, 0);
1162 for (
unsigned int k = 0; k <
data.n_shape_functions; ++k)
1164 const unsigned int comp_k =
1165 fe.system_to_component_index(k).first;
1166 if (fe_mask[comp_k])
1167 for (
unsigned int j = 0; j < dim; ++j)
1168 for (
unsigned int l = 0; l < dim; ++l)
1169 for (
unsigned int m = 0; m < dim; ++m)
1170 for (
unsigned int n = 0; n < dim; ++n)
1171 result[fe_to_real[comp_k]][j][l][m][n] +=
1172 (fourth[k][j][l][m][n] *
1173 data.local_dof_values[k]);
1177 for (
unsigned int i = 0; i < spacedim; ++i)
1178 for (
unsigned int j = 0; j < spacedim; ++j)
1179 for (
unsigned int l = 0; l < dim; ++l)
1180 for (
unsigned int m = 0; m < dim; ++m)
1181 for (
unsigned int n = 0; n < dim; ++n)
1183 tmp[i][j][l][m][n] = result[i][0][l][m][n] *
1184 data.covariant[point][j][0];
1185 for (
unsigned int jr = 1; jr < dim; ++jr)
1186 tmp[i][j][l][m][n] +=
1187 result[i][jr][l][m][n] *
1188 data.covariant[point][j][jr];
1192 for (
unsigned int i = 0; i < spacedim; ++i)
1193 for (
unsigned int j = 0; j < spacedim; ++j)
1194 for (
unsigned int l = 0; l < spacedim; ++l)
1195 for (
unsigned int m = 0; m < dim; ++m)
1196 for (
unsigned int n = 0; n < dim; ++n)
1198 jacobian_pushed_forward_3rd_derivatives
1199 [point][i][j][l][m][n] =
1200 tmp[i][j][0][m][n] *
1201 data.covariant[point][l][0];
1202 for (
unsigned int lr = 1; lr < dim; ++lr)
1203 jacobian_pushed_forward_3rd_derivatives[point][i]
1206 tmp[i][j][lr][m][n] *
1207 data.covariant[point][l][lr];
1211 for (
unsigned int i = 0; i < spacedim; ++i)
1212 for (
unsigned int j = 0; j < spacedim; ++j)
1213 for (
unsigned int l = 0; l < spacedim; ++l)
1214 for (
unsigned int m = 0; m < spacedim; ++m)
1215 for (
unsigned int n = 0; n < dim; ++n)
1217 tmp[i][j][l][m][n] =
1218 jacobian_pushed_forward_3rd_derivatives[point][i]
1221 data.covariant[point][m][0];
1222 for (
unsigned int mr = 1; mr < dim; ++mr)
1223 tmp[i][j][l][m][n] +=
1224 jacobian_pushed_forward_3rd_derivatives[point]
1227 data.covariant[point][m][mr];
1231 for (
unsigned int i = 0; i < spacedim; ++i)
1232 for (
unsigned int j = 0; j < spacedim; ++j)
1233 for (
unsigned int l = 0; l < spacedim; ++l)
1234 for (
unsigned int m = 0; m < spacedim; ++m)
1235 for (
unsigned int n = 0; n < spacedim; ++n)
1237 jacobian_pushed_forward_3rd_derivatives
1238 [point][i][j][l][m][n] =
1239 tmp[i][j][l][m][0] *
1240 data.covariant[point][n][0];
1241 for (
unsigned int nr = 1; nr < dim; ++nr)
1242 jacobian_pushed_forward_3rd_derivatives[point][i]
1245 tmp[i][j][l][m][nr] *
1246 data.covariant[point][n][nr];
1262 template <
int dim,
int spacedim,
typename VectorType>
1264 maybe_compute_face_data(
1265 const ::Mapping<dim, spacedim> &mapping,
1266 const typename ::Triangulation<dim, spacedim>::cell_iterator
1268 const unsigned int face_no,
1269 const unsigned int subface_no,
1271 const typename ::MappingFEField<dim, spacedim, VectorType>::
1280 const unsigned int n_q_points = output_data.boundary_forms.size();
1289 for (
unsigned int d = 0; d != dim - 1; ++d)
1291 Assert(face_no + cell->n_faces() * d <
1292 data.unit_tangentials.size(),
1295 data.aux[d].size() <=
1296 data.unit_tangentials[face_no + cell->n_faces() * d].size(),
1301 data.unit_tangentials[face_no + cell->n_faces() * d]),
1309 if (dim == spacedim)
1311 for (
unsigned int i = 0; i < n_q_points; ++i)
1319 output_data.boundary_forms[i][0] =
1320 (face_no == 0 ? -1 : +1);
1323 output_data.boundary_forms[i] =
1324 cross_product_2d(
data.aux[0][i]);
1327 output_data.boundary_forms[i] =
1328 cross_product_3d(
data.aux[0][i],
data.aux[1][i]);
1344 for (
unsigned int point = 0; point < n_q_points; ++point)
1349 output_data.boundary_forms[point] =
1350 data.contravariant[point].transpose()[0];
1351 output_data.boundary_forms[point] /=
1352 (face_no == 0 ? -1. : +1.) *
1353 output_data.boundary_forms[point].norm();
1362 cross_product_3d(DX_t[0], DX_t[1]);
1363 cell_normal /= cell_normal.
norm();
1367 output_data.boundary_forms[point] =
1368 cross_product_3d(
data.aux[0][point], cell_normal);
1374 for (
unsigned int i = 0; i < output_data.boundary_forms.size();
1379 output_data.JxW_values[i] =
1380 output_data.boundary_forms[i].norm() *
1381 data.quadrature_weights[i + data_set];
1386 const double area_ratio =
1388 cell->subface_case(face_no), subface_no);
1389 output_data.JxW_values[i] *= area_ratio;
1394 output_data.normal_vectors[i] =
1396 output_data.boundary_forms[i].norm());
1407 template <
int dim,
int spacedim,
typename VectorType>
1409 do_fill_fe_face_values(
1410 const ::Mapping<dim, spacedim> &mapping,
1411 const typename ::Triangulation<dim, spacedim>::cell_iterator
1413 const unsigned int face_no,
1414 const unsigned int subface_no,
1415 const typename ::QProjector<dim>::DataSetDescriptor data_set,
1416 const typename ::MappingFEField<dim, spacedim, VectorType>::
1420 const std::vector<unsigned int> &fe_to_real,
1424 maybe_compute_q_points<dim, spacedim, VectorType>(
1430 output_data.quadrature_points);
1432 maybe_update_Jacobians<dim, spacedim, VectorType>(
1433 data_set,
data, fe, fe_mask, fe_to_real);
1436 const unsigned int n_q_points =
data.contravariant.size();
1439 for (
unsigned int point = 0; point < n_q_points; ++point)
1440 output_data.jacobians[point] =
data.contravariant[point];
1443 for (
unsigned int point = 0; point < n_q_points; ++point)
1444 output_data.inverse_jacobians[point] =
1445 data.covariant[point].transpose();
1447 maybe_update_jacobian_grads<dim, spacedim, VectorType>(
1448 data_set,
data, fe, fe_mask, fe_to_real, output_data.jacobian_grads);
1450 maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
1456 output_data.jacobian_pushed_forward_grads);
1458 maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
1464 output_data.jacobian_2nd_derivatives);
1466 maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
1474 output_data.jacobian_pushed_forward_2nd_derivatives);
1476 maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
1482 output_data.jacobian_3rd_derivatives);
1484 maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
1492 output_data.jacobian_pushed_forward_3rd_derivatives);
1494 maybe_compute_face_data<dim, spacedim, VectorType>(
1495 mapping, cell, face_no, subface_no, data_set,
data, output_data);